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THE THEORY OF DUAL CONTROL, III 


A. A. Fel'dbaum 

Translated from Avtomatika i Telemekhanika, Vol, 22, No, 1, 
pp. 3-16, January, 1961 

Original article submitted April 18, 1960 


We shall study two examples of the synthesis of an optimal system for dual control, In the first 
example, we study a system of automatic stabilization with a linear object in which the per- 
formance criterion is the mean-square error. In the second example, we study an automatic sys- 
tem which searches for the minimal output value for an object that has a parabolic character- 


istic, The analysis is performed on the basis of the theory presented in the preceding papers [1,2] 
in this series, 


Remarks on the Synthesis of Optimal Dual-Control Systems 


In [1,2] we formulated the statement of the problem for the case of synthesizing optimal dual-control sys- 
tems, and presented a general method for solving this problem, In the analysis below we study examples of the 
application of the general theory to individual problems, 


The method of solution consists of formulating a series of functions y p~, (k= 0, 1,..., 0), and finding 
those values u,y- , of the control input which minimize the functions y p»-,. Assume the minimal value of the 
function Yp - , with respect to up - _ is denoted by y§-,. The function yp - -; is formed from y§- , by inte- 
grating the output value yp - , of the object which is measured by the control unit, and by adding the function 
Gn-k-1 Thus, the scheme for the basic computations involves a series of alternate integrations and minimiza- 
tions which is clearly depicted in Fig. 1.* The results of the computations are the optimal inputs uf . 


In order to find the optimal control input uf at the initial instant t = 0, it is necessary to pass through the 
entire chain of computations shown in Fig. 1. At the instant t = s the volume of the computations required for 
determining u$ is reduced, since it is necessary to perform only a portion of the chain of computations shown in 
Fig. 1— from the end of uf, to the determination of ug. The latter quantity is determined as a funcuon of the 
vectors U._ y, X§, and ys- 4; i.e., it is determined as a function of all the preceding values of the quantities at 
the input and output of the control section of the system, These values are memorized in the control section, 


Depending on the complexity of the computations performed in the chain depicted in Fig. 1, the following 
cases may be encountered; 


a, The entire chain of minimizations and averagings can be performed (either exactly or approximately) 
analytically, In that case, it is possible to compute the functions 


in advance and to introduce the ready control law in the form of (1) into the control section A of the system, 
This case is treated in the first example in this paper, However, such a possibility is encountered only in the 
simplest examples, 


* Figure 1 indicates integration with respect to xj, since, in the subsequent examples, hg = 0 and yj = xj. 


| 
| 


averaging 
(up) | 
(us) | 
miny, 


Fig. 1. 


b, It is impossible to perform the entire chain of computations in advance, and to obtain ready formulas 
of the type (1), This is the more general case, Such a problem is treated in the second example, In this case, 
the theory provides the possibility of achieving the synthesis of an optimal dual-control system, Here we design 
a control section A for the system, such that it automatically performs (either completely or partially) the com- 


putations shown in Fig. 1. 


The Synthesis of the Optimal Control Section for the 


Automatic Stabilization System 


The block diagram of the system is shown in Fig. 2, The output quantity x of the controlled object B is 
applied to the input of the control section A of the system, The stand- 
ard input x* (which is assumed constant), is applied to the other in- 
put of the section, The output quantity u from the section A of the 


g 


Fig. 2. 


system consists of the control input, which is applied to the input of 
the object B, together with a random constant bias ! and a random 
noise g. The values gg of the noise at the instants t = s are a series 
of independent random quantities with the same distribution density 
a(gs). The a priori distribution density P(u) = Pg) is specified by 
the random quantity 1. The object B has a linear characteristic of 
the following form: 


= -+- +p. 
The specific loss function is determined by the expression 


The over-all loss function is written as 


n 
Way W,= 
s=0 


n 


> (z,— 


§=0 


(2) 


(3) 


(4) 


where the number of cycles n is specified, It is necessary to synthesize a control section A, such that the risk R 
(the mathematical expectation of the function W) is minimal. 


In accordance with the general theory, it is required to find the probability density P(yj|#, i, uj). Since, 


in this problem, y = x, # is a scalar quantity, and the specified probability density is independent of i, we shall 
write it as P(x;| #, uj). It is obvious that, in this problem, 


P (x |p, ui) = 9 — ui — p). 


(5) 


In order to determine the optimal strategy, it is necessary to find the function «)(0 = k = n),which is 
given by the expression 


k 
ae = \ WP, (u)T] P wi) (xe, = 
i=0 


k 
i=0 
We shall study the case where Po(u) and q(g,) are normal distributions with zero values, Then 
average 
| 263) ’ q (8s) | (7) 
From (6) and (7) we find 
+00 
| \ (2" —2y}*exp | (8) 
20% 
We shall make use of the formula ([3]}, p. 188) 
attidg dq? (9) 
where a is a positive integer, and p and q are rational, 
For a = 1, we obtain yo aS 
14/8 
da = (1 + (10) 


We shall find the integral I, in the square brackets of formula (8): 


\¢ — exp {— 28 | da | exp{ zat a) 
where 
Ip — = 2, —u,—p=Ik. (12) 
Making use of formula (10), we find 
I, = V 2x [1 —u, —p)*) (13) 
from (11) and (12), 
Then 
(2m) 2 | i=0 


The first of the integrals in the square brackets of (14) can be found if we use the formula ({3], p.185) 
+00 


\ de = (15) 


The second integral in the square brackets can be determined by assuming # + u,— x* =z. Finally, 
the formula for c), is obtained in finite form. 


We shall assume 


Uy — = Wp, >) — 41) = 
i=0 
28’ Py , (16) 
1 
a, = , by 3a 
(2n) Ve, * (e,)? 
Then 
a = fag + by + (wpe + {— + a7) 


The quantity uk, which is contained in wy, is present in only one of the multipliers, and is not included in 
the exponential term. This fact greatly simplifies the minimization of a) with respect to u,. 


We shall now perform the series of computations indicated in Fig. 1. In formula (17) we shall write k = n 


and find the minimum of the function y , = %,, with respect to up. It is obvious that the condition governing 
the minimum is written as 


whence we find the optimal value uj for the control input: 
n—1 
— 4) 
(19) 
Then we find 
= (on + by) exp {— 4. (20) 
n n 20% dofe,, 
After that we find the function 7 p - ; (cf. [2]: 
+00 
Yn—1 = + \ din = 
= 
(21) 
+00 
(an + 6,) exp {— 268 + dans. 


+ = 0, (18) 
4 


We denote the integral in this formula by Jn-,. Using formula (15), we determine this integral: 


(22) 


x 


2% 235 235 + 


204 
This quantity is positive, since n = 1 and =z > 0. It should be noted that Jn is independent of up - y- There- 


fore, in formula (21), only the first term oy - , depends on up - 4; this facilitates minimization of yp -, with re- 
SPECE CO Up ge 


In this case, only - depends on up - and is independent of up - 4; this means that only 


the operational risk for t = n— 1 depends on up - ;, whereas the adaptive risk is independent of the quantity uy - , 
even though it is nonzero, Thus, the system under study is neutral (2), and in it the process of adapting to the 
object proceeds identically for any value of up - ;, and, in general, for us (0 = s < n), 


Minimizing ©, - , with respect to up - y, we find the condition governing the minimum in the same form 


as Eq, (18): 
Wn—1 + 233 = 0, (23) 
whence it follows that 
n—2 
— u;) 


Reasoning analogously, we find that, in general, y§ is independent of ws, and depends only on the differ- 
ence Xs-4— Us -y. Therefore, the substitution z = x,~4— us~, liquidates dependence on ug. , in the integral 
Js- ye Therefore, for any s— 1, only &,~ , depends on ug - 4, and the formula for determining u§ is written as 


n° = — (25) 


= 0, 3,..- a). (26) 


This formula has a simple meaning, If the noise g, were to be equal to zero, then (Fig. 2) we would need 
only two values (namely x; and wy) for determining the bias ; here # = xj— wy. It is obvious that, after this 
measurement, the value of uj, y, Which must be introduced, is equal to x*— pt = x* — (xj — uj). If the noise g, 
exists, then, for large s, it is natural to determine # in the form of the arithmetic mean of all the measurements 

s=1 
(x; — 4,) 
of uj), — Then we introduce the value 
Pay~ 


= (An + bn) x 
V 
| ide.) 
In this formula, we encounter the quantity o. 
1 1 1 
| 
n 
(=; — uj) | 
) | 
Substituting the value of €, from (16), we find the optimal control law in the following form: | 
s—1 | 
| 
g | 
s+ | 
| 


Tag" 


as) j4s Us = —pay= — (27) 
! few) The exact formula (26) differs from formula (27), 

5s 7 " ; which is derived from simple considerations, only for 
| smalls. For sufficiently large s , the two formulas coin- 
k cide, 
This example is the simplest one, and is illustrative 


in nature, since the measurement of the difference (xj — u4) 
Fig. 3. makes it possible to determine the value of # + g., and 
then to determine the bias as if it had been measured 
directly with the error g,. The same result would be obtained in the problem under study in the case of a non- 
linear, mutually single-valued characteristic for the object. However, even in linear systems, the synthesis of 
an optimal control section A is appreciably complicated when two coordinates 1, and [l, of the vector # are 
present, if the characteristic of the object has, for example, the form 


Ls = (Us py + (28) 
The general theory cited in [2] makes it possible to solve this problem in principle, 


From the optimal control law (26) it is not difficult to formulate the block diagram of the control system 
(cf, Fig. 3), The sign (—1) denotes an inverter (i.e., a sign-changing block), The letter T denotes a block that 
produces a delay of one cycle, and the letter C denotes an adder with storage. The result of the summing (it is 


assumed that n is not so great that the adder overflows) is the quantity }) (xi — ui), which is subjected to 
multiplication by [s + (2) ] in the block a(s), and then is subjected to sign inversion, 


It is of definite interest to extrapolate this scheme to the case of a continuous system, We shall assume 
that the duration of a cycle is At = t/s, and shall make s go to infinity while maintaining t = const, Then 
At + 0, Assume that we have the problem of synthesizing a control system for the case where the noise g(t) is 
stationary white noise with the spectral density Sy. However, white noise is not the limit of a series of independ- 
ent random quantities with finite dispersion, since its dispersion is infinitely great: 


2n 


x 
6 = \ Sedo = 


Therefore, it is possible to obtain white noise in the limit for At - 0 from a series of quantities g, only in 
the case where the dispersion of these quantities is correspondingly increased, Assume the correlation function 
for the noise is written as 


K,(t) = — = a(t — for —At<t< At; 


(29) 
K,(t) =0 for |t|> At. 
Then the spectral density is 
at 
S,(@) =2 \ Ke(t) cos = 
0 
202 
sin wAt — [wAt sin wAt -++- cos @At — 1}. (30) 
For At 0, the function tends toward ofAt: 
lim S,(@) = Tim of At = = lim At = Sp. (31) 


Ato At-»o 


14) 


9) 


0) 


31) 


4 
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7 
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Fig. 4. 
g 
Fig, 5. 


We shall rewrite formula (26) as follows: 


For At + 0 the quantities in this formula tend to- 
ward the following limits: 
s—1 
—> Sy, (a — wy) At —u) de. (33) 
i=o 9 


Since sAt = t, it follows that, in the limit, we obtain 
the following optimal control law from (32): 


t t 
2 


where a = S/o, The second term in this expression can 
be obtained at the output of a filter with the transient re- 
sponse 


i t 
Hout = = (35) 
0 

if we apply the quantity uj, = x— u to its input.* This 
transient response is not an exponential, The transfer func- 
tion K(p) for such a filter is expressed in terms of an in- 
tegral exponential function, and cannot be represented in 
the form of a finite sum of fractionally rational functions, 
However, it is possible approximately to replace the filter 
having the transient function (35) with an inertial section 
having the time constant T * 1,1 a, Then a system which 
is close to optimal will be of the form illustrated in Fig, 4a. 
Converting this system to an equivalent one, we successive- 


ly obtain the systems shown in Figs, 4b and 4c, The latter variant is a section which introduces the integral into 
the channel of the system, and has an inertial section in the feedback loop, It is of interest to compare these sys- 
tems with similar ones which were obtained on the basis of completely different concepts by I, B, Reswick [4] 


and O, T, M, Smith [5]. 


The Synthesis of the Optimal Control Section for an Automatic Search System 


We shall study an example of applying the theory of dual control for synthesizing an optimal automatic 
search system, The block diagram for the system is shown in Fig, 5, Assume that the characteristic for the ob- 


ject B is written as 


It is required that the value of x prove to be as close as possible to the minimum as a result of automatic 
search, The constant bias yp at the object's input is not known a priori, and is a random quantity with a speci- 


fied a priori probability density P{ 1)=P (pu). The noise gg.» just as in the preceding example, consists of a series 
of independent random quantities with a zero average value and the same probability density P(g.) = a(gs). 


*In determining the transient function, it was assumed that uj, = 1 = const, 


A s—1 
) >) — ae 
sAt + 
B 
) 
z 
t 

= 

1 


wy) 919s) In order to determine the required value of 


a b o the control input, we establish n trial values of us at 
% the instants s= Then, fors =n, 
we establish the final working value u = up, after 
Fig. 6, which search ends, It is required to establish the 


series of trial values u, (s = 0,1,...,m- 1), and 
the working value up in such a way that the mathematical expectation for the deviation of xp from the minimum 
is minimal, Such a search will be an optimal search performed according to the extrapolation method, 


The loss function for the problem under study will be written as 


W, = 0,2, = (us + +p)’, (37) 
where 
i, s=n. (38) 


The problem consists of determining an algorithm for the control section A which is such that the risk R = 
= Ry = M{Wp} will be minimal, Therefore, for any sth cycle (s < n) the entire risk is the adaptive risk, and the 
operational risk is equal to zero, Moreover, for s = n, the over-all risk is the operational risk, Further, the func- 
tion &, is equal to zero fors < n, The quantity 

n—1 
R= Ry, = (tn un) P x 
2 (zp, Uns Xn—)) i=o 


PidQ (an, Uns Xn—1) 


i=o (39) 
is to be minimized, 


In accordance with general theory, it is necessary to study the auxiliary function 


n—1 


(tn |p, tn) [] P (ail dQ (xn, (40) 
Q(Xp, i=0 


If the quantities and g, are normally distributed, then it follows that, from Eq. (40), it is possible to ob- 
tain the function &,y in finite form by integration, However, even in this case, only the first minimization (with 
respect to up) is not beset with difficulty, The subsequent integrations and minimizations already cannot be per- 
formed analytically, For other distribution laws, P(t) and q(g,), the situation is complicated even at the first 
stage of the computations, In the analysis below we shall study the case where the densities P() and q(g,) corre- 
spond to uniform distribution (Fig. 6), Assume also that the condition b < a is satisfied. The meaning of the 
parameters a and b is clear from Fig. 6, Assume that 


Ns = Us + (41) 
Then = and 2n dn, = From this it follows that* 
Pr (z’ << 2’ 4- dz,) = 


dz, | 7 dz, (42) 
= <n <V2 +5) +Pr(—Ve 
Pr(e’ 4-day) = Pr (V2 — 0, < VE — + + 


dx 
* The letters Pr are the abbreviation for “probability.” 


P (xj\p, ud) = V —ui—p) + — 
Assume 
Then 


Substituting this expression and the value of P(j) into (40), we find 


(44) 


(45) 


(46) 


(47) 


(48) 


Here we should take into account the fact that dx = 20 pdép. Formula (47) can be written in the follow- 
ing form: 
"7 
where 


(Ens +H) don. 


In the latter integral, the upper limit may be finite, since Q, = 0 fort, > a+2b, 
The computation of the integral I,(u) yields the formula 


Tn(p) 


In = + (tn 


Therefore, 
(I I ats uj + p) dp + 
+ + Qt dp). 

We shall assume that 


b n—1 n—l 


b n—1 


im 
Then formula (51) will become 


i 
1,7” 
n en 


(49) 


(50) 


(51) 


(52) 


(53) 


(54) 


Thus, 
im 
) 
- 
1 
— 
) 
9 


From object 


Q% 

| Bi -@ 

Sn Sn fou} 
MQ, | Q,--- [Q,_, 

Fig. 7. 
and the corresponding value of a5: 


The expressions Ip n, 1; ,n» 1z,n in the form of formulas are so complex that it is simpler to determine the 
integrals in the actual control section A, The further integrations and minimizations in formula (55) can also 
be performed only within A, since the formulas cannot be obtained in finite form, However, the over-all se- 
quence of computations can be realized by a computer device whose block diagram is shown in Fig. 7. This is 
the block diagram for an optimal control system, The blocks in this system can, in principle, be either analog 
or digital, and the diagram proper can be replaced by its equivalent program on a universal digital computer. 


In the block diagram shown in Fig. 7 we have used the following notation: B are rectifiers which are con- 
trolled by a control and synchronization block not shown in the figure; p are generators which produce “tri- 
angular® waves at various frequencies; V are the modulus blocks for which the output quantity is equal to the ab- 
solute value of the input quantity; S; are integrating blocks (a more detailed diagram of the integrating block is 
shown in Fig, 8); AO; are automatic optimizers with respect to one variable; II are the product blocks; DU is a 
dividing unit; BII, is the eg block for the uj values; BIlz is the memory block for ¢; values; ¥ is a block 
whose output quantity is €;= ¥x;, where xj is the input quantity, 


10 


5) 


The generator I'p,, applies a voltage to the bus  , which varies over the limits —b = » = b according to 
a "triangular" curve; here | = const. The blocks Qo, Qy, « « » Qgs « » Qn=4 form the corresponding func- 


tions in accordance with Eq, (45), As an example, Fig, 7 shows the block diagram for the section Q, in more de- 
tail; here Q, is the sum of the output signals from two rectifiers B, Each of these output signals is equal to $a 
or zero, depending on the conditions 


The terms q in expression (45) are obtained at the outputs of the rectifiers, mene 
All of the functions Q, (s = 0,1,...,0— 1) are applied to the block II, where the product | [ @, = 7Q, 
i=0 


is produced, The resulting product, and also the products #1IQ; and # *1Q;, are applied to integrating blocks 
where integration with respect to ! is performed during a half-cycle of the voltage from the generator I’ p,, . The 
resulting values Ip n, Iy,n,. and I, n are applied to the input of the block, which computes the function 
2 

u®To.n + Unlin + Tan (57) 
{i.e., that factor in the function ¢&p = 7, which depends on up in accordance with formula (53)), At the output 
of the dividing unit DU, the normal value of ug as computed according to Eq. (54) is determined; the quantity 
us is applied to the input of the block y,. Therefore, at its output, we obtain the value 73° corresponding to 
Up = ug. The function y, is integrated in the block Sp - , during one half-cycle of the generator T'pp - ;, i.e., 
it is integrated with respect to the coordinate €p-,, which, for the time being, arrives from the generator Ipp - ,. 
As a result, we obtain the function 7-4. In order to integrate we must first multiply a% [cf, (55)] by 
26n- den -4, and then integrate, Therefore, only the square bracket in formula (55) must participate in the in- 
tegration, The multipliers 1/¢j must not be taken into account, 


The output signal yj - 4 of the block Sp - , is transmitted to the input of the automatic optimizer AOp - 4, 
which chooses the value of uj ., in such a way as to ensure the minimum for yj. 3. The next integration over 
a half cycle of the generator I pp -» in the integrator Sp-, makes it possible to obtain the function yj -», etc. 
The chain of optimizers automatically selects the values of uf, uj, «+ , UA-4, Which minimize the functions 
Yor Vis ae Pespectively, 


The value of uf chosen in this manner is applied to the memory block BIly and stored in the cell uy. From 
this it is transmitted to the object at the instant t = 0, As soon as the value u, is memorized in the BIly and is 
applied to the object, the computation of up can be ceased, and from then on us is applied to the block Q, only 
from the memory block, Analogously, at the instant following the memorization of the value u,, the latter be- 
gins to be applied to the block Q, only from the memory block BIly, etc, Gradually, during the process of opera- 
tion of the network, all of the cells uj in the memory block are filled, 
Before the last cycle, the cell up is filled, In exactly the same way, 
all of the cells ¢j in the memory block BIiz are gradually filled, As 
soon as the value ¢; has appeared at the output of the object at the 
instant t = s, it is stored in the sth cell of the £, memory block, and 
from that instant the block Q, is subjected not to the output signal of 
the generator Ips, but to the constant value ¢, which is measured by 
the block, 


Figure 8 shows the block diagram of the integrator S. During the first half-cycle, the upper branch of the 
integrator I, operates (starting from zero) and accumulates the value of the integral. During the second half- 
cycle, the lower branch operates with an analogous integrator I,, and during that time the accumulated reading 
is taken from I, and stored in the output memory cell CII, Thus, we obtain the required value of the integral 


at the output with a lag equal to a half-cycle of the corresponding generator. 

The network in Fig. 7 is very complex, Its basic shortcoming, however, is not its complexity, but the 
necessity of an extremely rapid operation of a number of sections, The number of memory cells in the network 
is small; the number of blocks which perform the basic operations is also small, but the fundamental frequency 
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must be very great. Assume, for example, that n = 3, If the generator 'p, operates at a frequency of f = 5 Mc, 
then, during a half-cycle (i.e., during a period T = 107" sec), the integrators S, produce the values Ip,n, Iy,n+ Ie,ne 
The quantity 75 is applied to the integrator Sp - , with exactly the same lag, For integration with respect to 
Sn-4, it is required that no less than 10-20 values of yj}, correspond to a half-cycle of generator Ipp-dT py). 
Assuming that 20 half-cycles of Tpy must be used for one half-cycle of the generator I’p,, we obtain a fre- 
quency of 250 ke for the latter. Thus, during one of its half-cycles (i.e., during a period T, = 2 + 107° sec), the 
value y3 is produced, In order to compute uf, - ; = u$, it is necessary to pass over the entire scale —b = u, = b 
while testing not less than 10 to 20 positions on the scale, Therefore, the computation of u$ requires a time 
equal to Tf = 20T, = 40 + 107° sec, Reasoning analogously, we find that integration for obtaining 7} requires a 
time equal to T,; # 20Tf = 8 * 10°‘ sec, and for obtaining uf it requires a time equal to Tf * 20T, = 1.6 * 107? 
sec, Integration for the purpose of forming 7 § requires a time T, = 20Tf = 32 + 10 sec, and a time T,* = 

= 20T, = 20 * 0,32 = 6.4 sec is required to obtain uf, Thus, even the very high carrier frequency of 5 Mc (and 
at that for a relatively small value n = 3) cannot assure a short time for determining the first value up. Of 
course, the subsequent values uj will be determined more rapidly. If we were to determine the value of uf, not 
according to formula (54) but by minimization of y,, then the time for determining up would prove to be still 
greater, 


This example demonstrates that, as far as possible, we should attempt an approximate solution of the pro- 
blem involving the analytical determination of the values Uj since the computation of these values in a com- 
puter requires an extremely high speed of response, Modern digital machines cannot insure the required speed 
of response when complex problems are solved, This example also indicates the importance of papers in the 
field of designing ultrahigh-speed computers as far as the synthesis of optimal control systems is concerned, 


However, we should keep in mind the fact that, for the principle of computer operation cited above, the 
operating time is of the order of a", where a = const, Therefore, as n increases, the time increases extremely 
rapidly, If we perform the greater _ portion of the computations in advance, and feed already formulated func- 
tions into the computer, then the operation is accelerated, However, under these conditions, it is necessary to 
reinember functions of n variables, and thus the volume of the memory will be equal to a number of the order 
of b™, where b = const. Only the approximation of functional relationships in the form of functions of a definite 
_ type with a comparatively small number of unknown parameters, combined with an acceleration of the operation 
of computers and the development of approximate methods for solving problems, will create the practical possi- 
bility of designing optimal dual-control devices in cases which involve any appreciable complexity, 
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OPTIMUM PROCESSES IN SYSTEMS WITH DISTRIBUTED PARAMETERS 


A. G. Butkovskii 

Translated from Avtomatika i Telemekhanika, Vol, 22, No. 1, 
pp. 17-26, January, 1961 

Original article submitted May 23, 1960 


* This article is devoted to the problem of the optimum control of systems whose motion is general- 
ly described by nonlinear integral relationships that correlate the command inputs with the output 
variables of the system to be controlled, It is shown how the obtained results can be applied to 
the solution of problems in the optimum control of systems described by partial differential equa- 

tions in the case where the equation solution is expressed by an integral relationship, 


Important results in the theory of optimum processes have been obtained by using nonclassical variation 
methods in investigating lumped-parameter systems that are described by ordinary differential equations [1,2]. 
However, there is a large number of problems where the systems to be controlled have distributed parameters, 
Such systems pertain, for instance, to a large number of production-line industrial processes, in particular the 
heating of metals in through-passage furnaces, the drying of strip and friable materials, continuous etching and 
deposition of coatings, agglomeration, distillation, etc, The problems in the optimum control of systems with 
distributed parameters are formulated in [3], 


Statement of the Problem 


Many problems in the optimum control of systems with distributed parameters can be solved on the basis 
of the results obtained in considering the following optimum control problem, 


ion Let n coordinates of the vector Q = Q(t) = [Q,(t), . . « » Qn(t)) satisfy the equation 
ty 
Qi = Q(t) = Kilt, t, ae dt, 
where 
K,(t, (ii, 2,...,n) 
d ; are certain given functions of the arguments, The control functions uy, = u,(t) (k = 1, 2,...,1), are sectional- 
ly continuous with respect tot. The control vector u = u(t) = (u(t), . . . , up(t)] belongs to a closed region © at 
956). any instant of time t contained in the t) = t = t, interval, 
It is necessary to find such control u = u(t) € 2, ty = t = ty for which Q(t,) = Q, (Q, is the assigned vector), 
and the assigned functional 
ty 
Qo = | (2) 


assumes the minimum possible value. 


Optimum Regime Conditions 


The necessary optimum regime condition which the control vector u = u(t) (t) = t = t,) must satisfy can 
be formulated by the following maximum theorem, 


Theorem 1, Letu= ult) €Q (t) = t = ty) represent such control for which Q{t,) = Q, , where Q, is the 
assigned vector, For the optimum control u = u(t), and the corresponding trajectory Q = Q(t), there should exist 
such a constant nonzero vector ¢ = (Cp, Cy, « « » » Cp) that Cp = 0 and that, for almost all t, tp = t = t, the func- 
tion 


c, 
Me 
| 
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(3) 


+ cof (t, Q(t), u) + c; Ky (ty, t, w) 


of the variable u, u €Q, attains its maximum at the point u = u(t), 


In certain problems, the point Q(t,) must belong to a certain multiplicity set M with an arbitrary dimen- 
sion that does not exceed n- 1, This is the so-called problem with a mobile right-hand end. 


If the point Q(t;)€M is known, then, on the basis of the above theorem, optimum control, as applied 
to the problem with a mobile end, must satisfy the maximum condition, An additional condition determining 
the point Q(t,) in the multiplicity set M is the transversality condition [1], Let Q(t,) €M be the optimum tra- 
jectory end, We shall say that the transversality condition is satisfied if the (cy, ... Cp) vector is normal to 
the plane S, which is tangent to the multiplicity set M at the Q(t,) point. The necessary condition for the opti- 
mum regime in the mobile right-hand end problem can now be formulated in the following manner. 


Theorem 2, Let u= u(t)i€n (t) =t = ty), for which Q(t,) €M, For the optimum control. u = u(t), and 
the corresponding trajectory Q = Q(t) in the mobile right-hand end problem, it is necessary that there exist a 
vector c = (Co, Cy, « » « » Cp) different from zero, which will satisfy the conditions of theorem 1, and that, more- 
over, the transversality condition at the point Q(t,) be satisfied. 

The proofs of theorems 1 and 2, which are based on the application of “acicular” variations, and on the 
properties of convex multiplicity sets [1], are given in Appendix I. 


In certain problems, instead of the system (1), systems of the following form must be considered: 


Q = Qi(t) = Kutt, (i= 4,2,...,2). (4) 


In this, the condition Q(t,) = Q, also must be satisfied for t = t,, and the functional (2) must assume the 
‘minimum possible value, This problem represents a particular case of the previous problem if the Kj(t,T, u) 
function in system (1) satisfies the condition 

K,(t, t, u(t)) =0 (i =1,2,....”), (5) 


where 


In this case, theorems 1 and 2 are also valid and, because of (5), the function H will assume the following 
form: 


i =1 


i=} 
Application of the Obtained Results 


As an example of the utilization of the obtained results in solving problems in the optimum control of sys- 
tems with distributed parameters, we shall consider two problems, 


The first problem pertains to heat-exchange processes where it is necessary to secure, in the minimum 
time, a temperature distribution in a solid which would be close to the assigned distribution by means of an ex- 
ternal temperature field, or to minimize in a certain time deviations (in any sense) of the temperature distribu- 
tion from a certain assigned distribution. 
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Assume that the q = q(x,t) function describes the temperature distribution in a solid in dependence on the 
space coordinate x (0 = x = x4) and the time t (0 = t = t,), where x, and ty are certain given values, Further, 
we shall assume that the heating equations are given by 


The boundary conditions are 
oq 
(9) 
The initial condition is 
q(x, 0) = 0. (10) 


The optimum control problem consists in determining a control function u = u(t), restricted by the condi- 
tions N; = u(t) = N, (0 = t = t), for which the integral of the square of the solid temperature distribution de- 
viation q, = q(x,t,) from a certain assigned constant temperature qg assumes the least possible value at the in- 
stant of time ty, i.e., 


x, 
I= — (2, ty) = min. (11) 


As is known from the theory of partial differential equations [4,5,6], the function q = q(x,t) that satisfies 
the conditions (7}(10) can be represented in the following form: 


q= (2, i)= t) u(t) dt, (12) 


0 
where ¢ is a known function of its arguments (the Green function method), 


The minimization of functional (11) can be reduced (see Appendix II) to the minimization of the func- 
tional 


Qo= \ (Qi (t) (ty, u(t) de (13) 


under the following conditions: 


t 
aw 


0 0 
where R and T are certain known functions, and y is a number, 
In order to find the extremal control function u = u(t), we shall apply theorem 1, 
The H function has the following form: 


t 
H = ( u(t) (t,t) (t)dt+-cy ty)u (t)+-es. (15) 


From the transversality condition (theorem 2), we have c, = 0; if we take into account that cy = 0, the 
function H has a maximum if the following condition is satisfied: 


| 
) | 


u(t) = My) x 


x sign {rR (t,t) —§ [7 (t, t) + T(t, t)] u dx}. (16) 
0 


Condition (16) has the form of an integral equation for the u = u(t) function. 
For the v = v(t) function, which is determined by the relation 


(Na+) + + (17) 
this equation is equivalent to the equation 


v(t) = sign (a(t) b(t, (x) dr), (18) 


where a(t) and b(t,r) are certain known functions which are expressed in terms of the known functions in Eq.(16), 
while |v(t)] = 1 (0= t= 


We shall assume that the expression under the symbol “sign” in Eqs. (16) and (18) does not become identi- 
cally equal to zero, and that there exists a solution v = v(t) which has a finite number of discontinuities in the 
(0, t,) interval, In this case, the following method can be used for solving Eq. (18), We shall first check whether 
the v(t)= 1 function or the v(t) =—1 function is the solution of this equation for 0 =t=t,. if not, we shall as- 
sume that the solution v = v(t) changes its sign only once at the point 6, (0 <@, < t,), and we shall calculate the 
expressions @, and 2 under the symbol “sign” in Eq. (18) as functions of time t and the parameter 9,. In this, 
©, = &,(t, 6,) is calculated under the assumption that v = 1 in the first interval, and 4, = ,(t, 0,) is calculated 
under the assumption that v = —1 in the first interval, Then, if there exists a solution v = v(t) with a single sign 
reversal, at least one of the equations 6,(6,,0,) = 0, &,(0,,6,) = 0 must have a real solution, If no such solution 
exists, the solution should be successively sought, first in the class of sectionally constant functions (equal to 1 
_ or ~1) with two sign changes, then with three, etc., by successively calculating the 4, and #, functions in de- 
pendence on the parameters 0,, 0,,63, etc, Finally, at a certain mth stage, as a result of the assumptions made, 
there must exist a solution v = v(t) with sign changes at m points 61, 62, ...,9m, i, there must be a real so- 
lution of at least one of the two systems of equations 

M, (6;, 0,,..- =O (i=1,2,...,m), 


“4 (19) 
(6;,9;,... , = 0 (é = 41, 2,...,m) 


(0< 01 << O2 <<. << 

Here, the = &(t,,93, » » » » 9m) function denotes the expression under the symbol *sign® in Eq. (18) in 
dependence on time f (0 = t = t,) under the assumption that v(t) = 1 in the first interval for 0 =t=@,. The 
, = &(t,4;, ~~.» 9m) function denotes the same under the assumption that v = —1 in the first interval for 


Let us now consider the optimum control problem as applied to the process of heat exchange between a 


stationary and a mobile medium, Such processes take place, for instance, in heating billets in through-passage 
furnaces, 


Let q = q(x,y,t) describe the temperature distribution in an infinitely wide strip with the thickness x,,which 
moves with the velocity v = v(t) = 0 in the y axis positive direction, The x axis is perpendicular to the plate 
plane, and the x coordinate varies within the 0 = x = x, interval, The heating of the material through the plate 
upper surface takes place in the 0 = y = y, interval, 


The heating equations have the following form: 


at O<r<m, 0<y<y, 0<t<4). (20) 
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The boundary conditions are 
aq 


Oz = (21) 
=% 
q (2, 0, t)=0. (23) 
The initial condition is 
q(x, y, 0) = qo(z, y)- (24) 


Here, a and & are positive constant coefficients, and the thickness x, of the solid can be a function of 


t 
n=y— \¥(p)dp. 


In this case, the optimum control problem consists in determining a control function u = u(t), restricted by 
the N, = u(t) = N, condition, for which the integral of the square of the deviation of the solid average tempera- 
ture, with respect to its cross section 


\ t) da, 


from the assigned temperature qg, assumes the minimum possible value at the point y = y, in the time interval 
from 0 to ty, is, 


[= \ (yi, t) dt = min. (25) 
0 
By the nondegenerate substitution of variables 
t 


n=y— \v(p)dp, t=f, (26) 


Equation (20) can be reduced to the one-dimensional heat-conduction equation, and the q = q(x,y,t) function 
that satisfies conditions (20}(24) can be represented in the following form [4,5,6]: 


xy t 
(2, t, (by — v(p) dp) d& + 
0 


If we assume that t, satisfies the condition v(p)dp < y, then, by taking (27) into account, the aver- 
age with respect to the temperature cross section 5 = Gy.) will be given by 


t 
7=9(y, t)=Gly, n+ (ee, t)u(x) dr. (28) 


Here, G(y,t) is a known function representing the first term in (27), averaged with respect to x, and k(t,T) 
is the ¢(x,t,7) function, which is averaged with respect to x. 
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t 
By denoting Q, (¢) = Jk (¢, tule) de and substituting (28) in expression (25) for y = y;, the problem will 
0 
be reduced to the minimization of the functional 


Qo = (2) (29) 


(where ¥ (r) is a known function], under the condition that 
t t 


We shall now compose the H function, Considering (6), we have 


t 


+ col¥ + ek (tr, u(t) + er. 
According to the transversality condition, c, = 0, cy = 0, 
It is obvious from (31) that H has a maximum if 


By substituting in this expression the Q,(T) value given by (30), we obtain an integral equation with respect 
to the function u(t) to be determined, the solution of which will define the extremum, 


It is obvious from (32) that, in the case where ¥(T) = Q,(r), the first term of the H function becomes 
identically equal to zero, and that the determination of the optimum control function from the maximum prin- 
‘ciple becomes impossible. However, it is readily seen from (29) that the functional Qgto be minimized is identi- 
cally equal to zero, i,e., it assumes the minimum possible value, According to the terminology given in [2], 
this case corresponds to the so-called “particular trajectory section, The optimum control u = u(t) is, in this 
case, determined from the solution of the Volterra linear integral equation of the first kind: 


¥(t)= t) u(t) dt. (33) 


The methods used for solving equations of the type given by (33) have been sufficiently treated in the 
literature [7]. 


Thus, on the basis of the above reasoning, we can find the optimum control inputs for a number of systems 
with distributed parameters. 


APPENDIX I 
Proof of Theorem 1, Let u = ult) EN (t) = t = t,) represent optimum control for which 


ts 
Qo= \ F(t, Q(t), u(t) at 


assumes its minimum value. 
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Let us determine the control function u* = u*(t), which is obtained ‘by varying the control given by u = u(t). 
For this, we shall select the instants of time T;,T2,...,71 which satisfy the inequalities <T,<T,<... 
e271 <ty, and which represent continuity points of the control function u(t), We shall select arbitrary non- 
negative numbers 5t;,.. , and arbitrary points Vz, Vz , belonging to the control region 2, Con- 
sider the half-intervals | given by 


edt, (i =1, 2,..., 0. (35) 


Assume that € is so small that the half-intervals (35) never intersect each other, and that they are con- 
tained within the t) = t = t, section; then the u* = u*(t) control will be determined in the following manner: 


If t does not belong to any of the half-intervals 1,,..., Ij , then 
u® (t) = u(t); (36) 
ifték Gi=1,2,...,2), then 


u*(t) = (36") 
Then, the variated trajectory has the following form: 
= Q(t) + eAQ(t)+..., (37) 


l 


AQ, (x) = > [K; (t, ty vj) — (x, u Ot, 
j=1 


ty 
AQ) = > [K, (t,t). 0) — (x, ty, Ot de + (38) 


i=1 j=1 


l 
+. >) (FG, Qe), — F (x), BE). 


Essentially, the [AQp, AQ(t,)] vector depends on the choice of the Tj and vy points, and the 54 (i = 1,2,...,2) 
numbers, (In the case where the same points Tj correspond to different variations, control variations similar to 
those used in paragraph 9 of [1] can be introduced.) Moreover, it can be considered that the same number of Tj 
and vj points have been taken for different variations, so that the addition of new Tj and vj points for which 
5t; = 0 would not change the control u*(t) that is being varied, 


Let us now consider the multiplicity K of vectors [AQy, AQg(t,)] with the origin at the point [Q),Q(t,)) in 
space (Qy, Q) which are obtained from the multiplicity of variations of the indicated type. We shall refer to the 
multiplicity set K by the term “cone of attainability,” It can be readily seen that K is a convex cone, Actual- 
ly, if the points s* and s* belong to cone X, the point s = \ 48° + A,s", where Ay =0, A, = 0, also belongs to coneK, 
since the point s is obtained by variation which is determined by the joining of T{ and vj for s*, and of r* and 

vf for s*, while 5t; for the point s is equal to 51; = Ay5t{ for Tf, and to tj =AgStf for rf. Since Eqs, (38) 

are linear with respect to tj, we find that the point s €K, Further, following a line of reasoning similar to that 
used in proving lemma 3 in [1], we find that the inside of cone K must not contain a straight line L with 
the origin at the apex of cone K, and which follows the negative direction of the Q, axis in the vector space 

(Qo, Q). Hence, it follows that there are such numbers C9, Cy, . + - » Cp for which the entire cone K lies in the 


n n 
- half-space b> 42, <0, and the straight line L lies in the half-space >) C_Q_ > 0, ise, the vector with the co- 


n a=0 
ordinates Cy = —1, = 0, » Cp = O lies in the lalf-space ¢,2, > 0. Consequently, cy = 0, whence we 
find that we 


a=0 


Let us consider the [AQy, AQ(t,)] vector, which is obtained as a result of varying the control function u(t) 
at the single point T, (2 = 1) and 5t,=1, By taking (38) into account, we have 
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n 
i, 


n 
F (41, Q(t1), —F (41, + > (ti. ta — Ky (ta, ta, <0. 
i=1 


Since this inequality is valid for any point vy €2, the function 


= Co 2 oF aq, * q(t. t, u) dt+ coF (t, Q(t), (41, l, u) 


must have a maximum for u €2, 


Proof of Theorem 2, Let us draw through the point (Qo, Q(t,)] in space (Qy, Q) a plane P that is 
parallel to space Q, We shall also draw a plane R with the dimension n (parallel to the Q, axis) through the 
plane S (with the n— 1 dimension in space Q)~which is tangent to the multiplicity set M at the point Q(t,) €M- 
and the point [Q,, Q(t,)] in space (Qy, Q), Consider the half-plane T of the plane R that lies in that half of the 
space divided by the plane P where the negative "end" of the Qy axis is located, Due to the fact that the Q = 
= Qt) (t) = t = t,) trajectory is the optimum trajectory, the attainability cone K must not contain any straight 
line L belonging to the half-plane T and having the origin at the point (Qo, Q(t,)]. In the opposite case, as at 
the beginning of the proof of lemma 10 from [1], we would arrive at a contradiction to the fact that Q = C(t) is 
not optimum, i,e., we could reach the multiplicity set M for a smaller value of Q). Cone K and the half-plane 


sect with cone T, Therefore, there must exis a plane given by the equation =0 (¢o<0), for which 

cone T lies in the half-space s 42, >0, and cone K lies in the half-space > ¢,2, <0. For instance, the 
a=0 a=0 

plane R, containing the half-plane T which, in turn, by definition, contains the plane S tangent to the multi- 

plicity set M, is such a plane, This proves that the transversality conditions are satisfied, 


APPENDIX II 
By substituting (12) in (11), squaring, and by expanding into a sum of integrals, we have 


Let us reverse the order of integration in the second term of the sum, and let us denote R(t, t) = 


@(z,t, t)dz andy = 2q,. Let us represent the square of the integral with respect to time in the third 


as a double integral, and also reverse the order of integration while denoting 7 (t, t) = 


\ (2, t) t) dx. Then expression (42) will assume the following form: 
0 


t t, t 

1 R (ty, t) u(t) dt + T (t, t) u(t) u(t) dt dt. 
0 

It is obvious that (43) can be written thus: 


¢ 
= 293+ T (t;, T) u(t) dt —yR (ty, dt. 


0 
corm 
‘ 
0 0 


2) 


3) 


A) 


By denoting Q: (¢) = \ T (t, t) u(x) dt and neglecting the constant term x,q%, we find that the minimiza- 


0 
tion of functional I is equivalent to the minimization of the functional 


ty 
= \ [Qu (t) — u(x) de 


Qo 


under the condition that 


t, t 
(er (46) 
0 
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CERTAIN PROBLEMS IN DESIGNING SYSTEMS WITH 
EXTREMAL SELF-ADJUSTING CORRECTIVE DEVICES 


V. N. Varygin 
Translated from Avtomatika i Telemekhanika, Vol, 22, No, 1, 
pp. 27-36, January, 1961 

Original article submitted July 6, 1960 


This article presents the derivation of equations of the extremum self-adjustment of corrective 
devices in linear systems for scanning oscillation frequencies that are comparable to the natural 
frequency of the basic control loop, These equations make it possible to investigate self-adjust- 
ing systems that are faster than those which can be investigated under quasi-steady~state condi- 
tions, The optimum relationships between the gain factors in self-adjustment loops are derived, 
Some of the results obtained in experimental investigation of such systems are given, 


If elements whose parameters change at random (the system to be controlled, the regulator, etc,) are pres- 
ent in the automatic control loop, it often becomes necessary to solve the problem of the system dynamic char- 
acteristic stabilization by suitably varying the parameters of the correcting devices (CD), The principle of CD 
extremal adjustment can be used for this purpose [1]. 


The schematic diagram of such a system is shown in Fig. 1. The input signal @ is fed to the basic system 
Wp and to the reference system W,, which is synthesized on the basis of the desired control performance criterion 


| 


and the available information on the input signal. The 
difference € between the reactions of the basic and the 
reference systems serves as the index of the deviation of 
the W), system characteristics from the assigned character- 
istics. The quantity €* can be conveniently used as the 
measure of this deviation, 


The CD adjustment parameters X;, X3,..., Xn 
are subjected to random or regular scanning oscillations 
with a rather small amplitude, We shall consider only 
regular sinusoidal scanning oscillations* [7]. 


Synchronous detectors SD,which are connected to 
the error signal quadrator output, separate the components 
which are synchronous with the scanning oscillations, and 
which modify the average values of the adjustment param- 
eters X;, Xg,... . Xp through the integrating elements 
until the system arrives at the state corresponding to the minimum é’, 


If the variations of the Wp system parameters are completely compensated by variations of the CD adjust- 
ment parameters, then, under steady-state conditions, the basic and the reference systems will be identical, 
Under transient conditions, the CD parameters will follow the variations of the basic system parameters with a 
certain dynamical error, 


Equations of the Self-Adjustment Processes 


The general first- approximation integral-differential equations describing the dynamics of the self-adjust- 
ment processes have been derived in [3]. These equations cannot be solved for the general case, and they are 
too complex for engineering purposes, Moreover, the above article considers the quasi-steady-state , self- adjust- 
ment conditions, which correspond to low scanning and self-adjustment rates in comparison with the processes 
taking place in the basic control loop. For this case, very simple equations, which make it possible to reduce 
investigations of the dynamics of such systems to the investigation of linear automatic control systems by ordi- 
nary methods, have been derived. 


The scanning rate must be increased if the self-adjustment system operating speed is to be increased, This 
would make it possible to broaden the transmission band of filters at the output of the SD multiplying elements, 
which basically determine the speed of self-adjustment processes, The present paper is concerned with the above- 
described system for the case where the scanning oscillation frequencies are comparable to the natural frequency 


For a small scanning oscillation amplitude and sufficiently slow changes in the values of the CD adjust- 
ment parameters, the system operator W = W,— Wp can be represented by the following fast-converging series 
(see Appendix): 


Fig. 1. 


W (D, Xo, Iny t) = W,(D, Xo, t) + 
> Lip COS ; (D, Xi, Xn, t) + 


ie] 


+ sin (D, Xn. t) + (1) 
i= 


P(D, Xn; t), 
where D is the differentiation symbol, w4 is the scanning oscillation frequency with respect to the ith parameter, 
W, and Wj are the system operator components which take into account, in the first approximation, the presence 
of scanning oscillations with respect to the ith parameter (see Appendix), and 


* The proposed investigation method can readily be extended to the case of periodic scanning oscillations of any 
form if they can be represented by a convergent Fourier series, 


= 


W,(D, Xi, Xo, Xa, t) W,(D) = W,(D, Xi, Xs, Xa; t). 


(The dependence on time is expressed by the fact that the X,, Xz,... , Xp values are considered as “frozen” at 
the instant of time t .) 


The P(D,X,,X2,. . . , Xn,t) operator includes the remaining terms of the series; since they include x4. fact- 
ors of higher power, i.e., terms which are small quantities of higher orders in comparison with Wj and W], they 
can be neglected, 


The system weighting function can be written in the following form: . 
K(X, X3,.., Xn, t, t) Xn, t, + 


+ >} cos wit Ki (X;,Xq,..., Xn, t,t) + (2) 
i=1 


i=1 


where KX eee ,Xp.t.7), Kj (X, Xe. eeer Xn t, Tt), and Xe, t, T) correspond to the 
operators W,(D, Xe, Xn, t), W,O, Xe. Xn» t), and wi(D, Xy Xn» t). 


Considering that the Ky, Kj, and K{ weighting functions are analytical functions with respect to the adjust- 
ment parameters, they can be expanded into series which have the following forms: 


j=1 


(OK 
Ky t) = Kg (sy!) + iis (3) 
j= 


OK; 
Ki (Xi, t,t) = Kj, (ty) + 2 + 
1 


where AX; = Xj — Xj, and Xjy is the adjustment parameter value corresponding to the steady-state conditions in 
the self- adjustment loops; the index denotes the substitution of the 


Xj = (j =1,2,...,m) 


values in the expressions for the weighting functions and their derivatives, 


According to design, the system must maintain the AX; (j = 1, 2, .. ., m) values close to zero; therefore, 
we shall hereafter consider only the first two terms of series (3), 


By substituting (3) in (2), and neglecting terms which are small quantities of higher orders, we obtain 


n aK n 
K (t, *) = Kg (t,t AX cos @t (t, t) + 
(t,t) or ( )+ 2 (az'), i+ 2 ar (t,t) te 


+ Ligh; Sin wt K;, (t,t). 
i=1 


If we assume that the variations of the average values of the adjustment parameters are inconsiderable 
within the range of large values of the basic loop weighting function K(t) [6], we obtain 


e(t)= | K (5) 


| 
n 
|_| 
r, 
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Considering that the input signal is a random function of time, we shall write the 0(T,)@(r,) product in 
the following form [4]: 


where R(T, T,) is the input signal correlation function, and f (T;, T,) is the centered random function, 
By substituting (4) in (6), and by taking (7) into account, we obtain 


=23 Lig COS wit S; 123 Tig Sin + R, (8) 
i=1 
where 
Si = j {Kor (t, + | \ Kr (t, ta) R(t, (9) 


while R comprises the constant component, different combination frequencies, and the centered random functions. 


The scanning oscillation frequencies must be chosen in such a manner that their combinations in the quad- 
rator signal do not produce a constant component in multiplication by the SD reference voltages, At the output 
of the ith SD multiplying element, we have 


u, = 00s w,fe* (t) = + Q,, (11) 


where Q; contains different, newly formed combination frequencies and centered random components,* 


We shall assume that the linear portion of the self-adjustment circuits completely filters out the Qi com- 
ponents, Actually, SD are resonance filters and, therefore, this condition can readily be satisfied by selecting 
the necessary operator Ws if the input signal @ spectrum does not have peaks near the scanning oscillation fre- 
quencies, Then, according to the diagram shown in Fig. 1, 


Xi = —kiW, (D) Si (i = 1,2... ., 00), (12) 
where 
= kiur?,. 
Let us write Sj in the following form: 
Si = + AX Sy, (13) 
j=1 
where 
Soi = \ \ Kor (t, Kir (t, Ta) R dt, (14) 
Sij = \ \ EE | Ki, (t, t2) dt, dt, 


Then the system of equations (12) in expanded form can be written thus: 


n 


DX, = — (D) [Sor + AX;S,5| 
1 


* Only one term of the €* series, namely the 2xjpcos witSj term, produces in multiplication by xj9cos wjt a con- 
stant component that can be separated by the Wf filter. 


(D)|Su+¥ 
1 


Under steady-state conditions for the self-adjustment system, all AX; (j = 1, 2,...,m) and all DX; (i = 
=1,2,..., 1m) will be equal to zero, Then, the adjustment parameter values corresponding to the extremum 
can be determined from the system of equations 


So =O (i= 14, 2,..., (17) 


By taking into account (17) and the equation DX; =. DXj; + DAX{, the characteristic equation of system 
(16) can be written in the following form: 


[Su | Sie Sin 
[ Se + | Sin (18) 


The forms of Eqs. (17) and (18) are similar to those of the equations obtained in [3] for the quasi-steady- 
state conditions. The equations become identical if we assume that the scanning oscillation frequencies are in- 
finitely small in comparison with the system natural frequency. In this case, 


9Ko(t, 1 de® 1 
The characteristic equation (18) can be written in the following form: 
+ kySy0 kySin 
+]... keSon = 0, (19) 
knSne knSan + I 
where 
(20) 


If there is an extremum (a minimum), all Ij values are negative [5], After determining the i (i = 
=1,2,...,) valves from the solution of the equation 


(21) 
the roots of the characteristic equation (18) can be found. 


The solution of Eq. (21) can be reduced to the investigation of the stability of the hypothetical system 
shown in Fig, 2, 


The characteristic equation of this system is identical to Eq, (21): 
—1,W; +1 =0. 
Thus, the stability condition for the system under considerationcan be defined (for lj < 0 for all i) as the 


stability condition for an isolated channel (Fig. 2) for lj = Imax. Actually, if the real parts of the roots of 
Eq. (21) are negative for Ij = Imax. they will be certainly negative for all ly < Imax. 


| 
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The control time will be defined by 


i 
Fig. 2. where A min is the root with the smallest absolute value which is found from 


Eq. (21) for lj = Imin (if the roots are real; if the roots are complex, the 
smallest value of the root real part must be used instead of min). 


If the difference between Imax and Imin is large, the processes will be delayed, since the Wf filter gain, 
which is determined from the stability condition for Eq. (21) for Ij = Imax, will yield insufficiently large root 
values in substituting = Ip in. By determining the optimum relationship between kp, » Kp, the ratio 
‘Imax/Imin can be reduced to the minimum,* 

Let us find the optimum relationships between the k, and k, coefficients from the condition for the 


Imax/Imin minimum for a system with two adjustable parameters, The characteristic equation for this case will 
have the following form: 


kySy +] 


= 0. 
+ 


It can be shown that the minimum of I,/1, with respect to k,/ kg will be obtained for ky/k, = Sg2/Sy, 
for kySyy = 
This is apparently the optimum relation also for a large number of adjustable parameters, i.e., for any n, 


= = kenSan (23) 
must hold, 


The performed experiments indicate that the use of condition (23) leads to a tenfold, and sometimes 
hundredfold, increase in the self- adjustment process speed in comparison with the case where the channel gains 
are equal, 


If the I43x and 1,4) values cannot be brought closer together even for the optimum kj value, the method 
of coupled control is recommended [1], 


Introduction of Scanning Oscillations into the Reference System 


In a number of cases, the introduction of scanning oscillations into the basic system either is inadmissible, 
or it entails considerable technical difficulties, Then, the application of scanning oscillations to the correspond- 
ing reference system parameter can be recommended, Such scanning has a number of important advantages: 


1, The scanning oscillation amplitude can be made much larger than the amplitude in the previous scan- 
ning variant, which would improve the self-adjusting system stability and noise-stability. 

2. No distortions due to scanning oscillations would be present in the main system during the signal trans- 
mission through the basic loop, 

3. By arranging the tracking of parameter changes in the basic system by the parameters of the reference 
system (the basic system model), and by introducing scanning oscillations into the reference model, the solutions 
of a number of purely extremal problems, which cannot be solved by ordinary methods, can be found, 


A disadvantage of the method of introducing scanning oscillations into the reference system is the neces- 
sity of designing reference systems with variable parameters. 


*In the limiting case, it can be reduced to unity. 


treg, sec In conclusion, it should be noted that, on the basis of the 
Xin. Obtained experience gained in practical work, the equations derived 
: from (22) experimentally under the assumption that the self-adjustment processes are 
3 45 45 slow in comparison with the basic loop operating speed are 
4 23 22 also valid with sufficient accuracy in the case where the basic 
5 13 12 loop operating speed is comparable to that of the self-adjust- 
ment loop, 
The author hereby expresses his gratitude to G. Yu. Polkanov for his great help in performing the experi- 
ments, 
Example, Let us consider the system for the transmission band stabilization in a very simple servosys- 
tem whose equation is given by 
ky = ke. (24) 


The gain is subjected to periodic scanning oscillations about a certain average value k = k,, 1.e., 


k = ka+ ky cO8 Wot (25) 
for ky << Ka. 


With sufficient accuracy, Eq. (24) can be written as j(Tg— Ty cos wet) + y = x, where 


he 
In this case, we shall use the following inertial element with the required time constant T, as the refer- 
ence system: 
W,(?) = (27) 
Wot D @ sin wt D (28) 
and 
ATD 
W (D) = W; (D) — Wa (D) = pie 
D (29) 
Considering (4), we can write 
t 
Ky FeV cos (Wet + @), (30) 
4 
= 
Hence, Sp, and Sy can be determined by using Eqs. (14) and (15), For an input signal @ = X;_, cos wt, 
we obtain 
(Ta — Tr 1 (1 — Trop) (1 — + Trop) + (31) 
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(1 — (1 — + + 47707) w* 
2 [AT + (1 — + (7 + 1) 
The value T,,, = T;, found from the solution of the equation So; = 0, has been substituted in expression (32), 
In this case, the self-adjustment process characteristic equation has the following form: 


(32) 


Su = 


33 

Su— (33 ) 

By substituting the expression for W¢ in Eq. (33), we can find all roots of the self-adjustment process char- 
acteristic equation, 


The oscillograms in Figs, 3a,b, and c illustrate the above example for the following data: 


1 sec; Tg= 0,1 sec; 


-1 1, 


W 1,26 sec™; w= 6,28 sec; 


1 
Xin = 3 v (for Fig. 3a); Xjp = 4 v (for Fig. 3b); Xin = 5 
(for Fig. 3c), 


The regulating time treg for the initial deviation from the extremum, calculated by using formula(22), and 
obtained experimentally, is given in the table, 


The possibility of adjusting two parameters is illustrated in the oscillogram given in Fig. 3d, which shows 
the transient processes in regulating the damping decrement and the natural frequency of an oscillatory element 
with the transfer function 


D? + 20 ED + 
for 
"t= Dit + 
and 


0.7, 1 sec! 


in the case of initial deviation of w and & g from war and Er. 


wel 
| Fig. 3. 


APPENDIX 


Determination of the Transfer Locus of a System 


with Extremal Self-Adjustment of the Corrective Devices 
The system equation with variable coefficients has the following form: 


(t) D® +... (t) D + ao (t)] y (t) = (t) D™ 4... + (t) D bo (t)) 2 (2) (34) 


L(D, t) y(t) = K(D, t)z(t); (35) 


where x(t) and y(t) are the input and output signals, respectively; a(t) and b(t) are certain known time functions, 
and D is the differentiation symbol, 


In [6], frequency methods, which are widely used in investigating systems with constant parameters, have 
been applied to systems with variable parameters, The H(jw,t) function, which is here used as the frequency 
response analog, is found by solving the following equation: 


(t)D™ +... + (t) D + (t)] H (jo, t) = K (jo, ¢), (36) 


where (1). and KGw,t), L(w,t) are obtained by substituting jw for D in the expressions for 
(j@)* 
K(D, t) and L(D, t). 


By using the method proposed in [6], the solution of Eq, (36) is obtained as a series with the following form: 


H (ja, t) = Hy (jo, t) + He(jo, t)+..., (37) 


which is a fast-converging series if the coefficients of Eq. (36) oscillate about their average values, which actual- 
ly takes place in the system considered in the main portion of the present article, 


For determining the terms of series (37), Eq. (36) should be rewritten thus: 


+... + a,D + ao] H = K + P(A), (38) 
where Oy is the average value of &,,(t) [or a constant approximating ©, (t)], and P(H) is the totality of variable 
terms in the equation, i.e., 

P (H) = H (a, —a,, (t)] D® +... + [a1 — (¢)] D + — a (#))). (39) 
The steady-state solution of the following equation with constant coefficients will represent the first ap- 
proximation of H(jw, th 
[a,D"+...+ + ao] Hy = K. (40) 
The H,, Hs, .. . values are obtained from the steady-state solutions of the equations 


[a,D" +... H, = P(H,_,) (41) 


In our case, if we consider the accepted assumption that the average values of the adjustable parameters 
change little in the region of large values of the basic control loop weighting function, and that the scanning 
oscillation amplitude is small, it is sufficient to determine only the first two terms of series (37), 


As an illustration of this, consider the above example, where 


K (jo, t) =1, L (jo, t) = jo (Tg — T, cos + 1. 
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For this case, Eq. (38) will assume the following form: 


Ta 1 To (i + D) 


Then, according to (40) and (41), ; 


(Tajo + 1° + + + Teas) (To +1) 


The next term of the series will be found by solving the equation 


Hy = T cos Wot 


T 
(T,D + 1) Hs = cos wol oF, 
and, since it will contain the factor Tj, it can be neglected. 


Every adjustment parameter in Eq, (1) corresponds to a pair of xjgcos wytWj(D, +,» Xp. t) and 
Sin wotWi(D, Xy, Xp, t) terms, the origin of which is clear from the above example, where 


jo w’ j@Ta 
(Taio +1)*+ [(Taj@ + 1)* + +1) 


The substitution of D for jc in (1) is possible, since W, and W/ represent the transfer functions of a system 
with constant parameters, 


t5=T,, W,= 
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DETERMINATION OF THE OPTIMUM VARIATION MODE 
OF THE USEFUL SIGNAL AND NOISE CARRIER FREQUENCIES 
IN DETECTION PROBLEMS BASED ON THE THEORY OF GAMES 


M. Yu. Gadzhiev 

Translated from Avtomatika i Telemekhanika, Vol. 22,No. 1, 
pp. 37-46, January, 1961 

Original article submitted May 16, 1960 


The conflicting situation arising when the useful signal is received against the noise background, 
in the case where the useful signal carrier frequencies and the noise can change within the as- 
signed frequency range, isconsidered here, The possible payoff function for the game model of 
the situation is determined, and the physical meaning of the payoff function is defined, The 
payoff function is approximated by expressions which make it possible to find the optimal stra- 
tegy and the game value by means of the existing methods, Concrete examples are given, 


The natural course to be followed in noise interference control is that of attempting to tune it out bychang- 
ing the useful signal carrier frequency. If, after the useful signal frequency has been changed, the noise carrier 
frequency is changed in such a manner that the most effective interference by noise is secured, a conflicting 
situation arises, where noise and the useful signal have opposite “interests.” 


Under certain simplifying assumptions, this situation can be represented in terms of the theory of games 
[1,2] by assigning the role of one (winning) player to noise, and the role of the other (losing) player to the useful 
signal, In the game corresponding to the conflicting situation under consideration, either of the players has only 
one personal move: the independent and simultaneous selection of his carrier frequency value from the closed 
interval [W,, W2]. 


Let us introduce the basic concepts and terms of the theory of games as applied to the problem under con- 
sideration: f, is the pure strategy of the first player (noise), which consists in the selection of the carrier fre- 
quency fy; f, is the pure strategy of the other player (useful signal), which consists in the selection of the fre- 
quency f2; ®(f3, f2) is the payoff function; € (f,) is the mixed strategy of the first player; n (f,) is the mixed 
strategy of the second player; € o( f ;) is the optimal mixed strategy of the first player; 1 9( f,) is the optimal mixed 
strategy of the second player; v is the game value, 


In operation, the useful signal receiver is approximated by a linear weighting function. In the general 
case, the noise state is assumed to be unsteady, 


Selection of the Payoff Function 
Let the useful signal g(t) become mixed with noise X(t), thereby forming the signal 


Z(t)= Xo(t) +o, (1) 
which arrives at the receiver, 
Let us denote 


X, (t) = X (t) exp i2nj,t = N (t) expi2n[f,t + a], (2) 
where N(t) and @ are the envelope and the phase of high-frequency oscillation (2), respectively; 


Ro (41, ts) =R (4, exp (4; ts) (3) 
is the noise complex correlation function; 


| 
| 


Go(t) = = +B) (4) 


is the weighting function of the useful signal linear receiver; h(t) and 6 are the envelope and the phase of the 
weighting function (4), respectively. 


We shall use the following expression for the payoff function: 


2 
O(/, fs) = =| Xo = 
s—-T 
= \ Ro (ty, ta) (tx) So (ta) dt (5) 
s—T s--T 
where (s— T, s) is the interval of integration by the receiver of signal (1); the asterisk denotes the complex con- 
jugate of the corresponding quantity, and the bar on top denotes statistical averaging. 


By taking into account (2) and (4), the payoff function (5) can be written in the following form® : 


/2) = \ \ R (ty, ta) h(t) h (ta) exp (t, — dt, dt, = 


—T 
ise 


8-—T s—T 
where Af = fz. 


Let us determine the possible physical meaning of the payoff function (5), We shall demonstrate that the 
payoff function (5) is a generalization of the case of unsteady-state noise, where the possibility of changes in the 
carrier frequencies to be effected by the players, the so-called “noise function,” is taken into account. This 
-function has been investigated by Woodward [3] for white noise and for the coincidence of the carrier frequencies, 


Let us consider the problem [3] of determining the “essential” parameter U of the g(t, U) signal. It is 
known [4] that, for extracting the useful information on U that is contained in Z, the optimal receiver must per- 
form the following operation: 


n(U,U’)=Re Z(t, U)dt = 
s—T 


= Re @ (t, U) g* (t, U’) dt + Re \ Xo(t) U’)dt, (7) 


s—T s-T 


where U® is the assumed value of the “essential” parameter U. The optimum weighting function g,(t) will be 
determined by solving the equation 


Fo (tas ta) (te, dts = @ (ta, (s) 


s—T 


*If the averaging is performed with respect to the phases of high-frequency oscillations (2) and (4), it can readily 
be shown that 


[Re (fi — = [Imo (i — fol? =F 


where Re and Im are the real and the imaginary parts of the complex quantity. 


lily 


while the solution is such that 


gn(t) op for s—T<t<s, 
\=0 for s—T>t>s, (9) 


Let us consider the terms in Eq. (7), where Re \ @(t, U) gt (t, U')dt= f (U, U’) ts the so-called signal 
8 s—T 
function, and Re \ Xq(t) (t,U’)dt= q (U’) is the noise function, 
s—T 


The U parameter is evaluated with respect to the signal function maximum (3), which is attained for U=U": 


max f(U,U')=m(U)= \ 9(t,U) U) dt. (10) 


s—T 


The error of the parameter U determination can be estimated with respect to the value of the noise func- 


tion mean square: 


= \ \ Ro (ty, ta) (ta, U’) Bn (toy U") dt, dt, = 


s—T s—T 


= 2[(Re = 2 (11) 
Expression (11) coincides with the payoff function (5)? 


Let us define more accurately the payoff function (6) for the particular case of steady-state noise by chang- 
ing the integration limits to (—oo, oo) in (5) in correspondence with formula (9): 


(Af) = \ R (a) H (a) {cos 2nAfa} da. (12) 


H (a) = h(t)h(t +) dt = K (f) exp (i2nfa)df = 
on \ | P exp df. 
if S(f)= exp (—i2nfa) da ts the spectral dent, then, with the new notation, 


Eq. (12) will assume the following form: 


(Af) = K (fy SN ay, (14) 


After the substitution of variables, Eq. (14) can be written as 


* By using (8) and (9), it can be shown that the signal-function maximum is equal to the mean square of the noise 
function, i.e., 


n (U') = m( 0). 


Moreover, the noise function mean square is equal to the generalized signal-to-noise ratio (5) in using 
the optimum filter, which is determined from Eq. (8). 


4 
ies, 
Here, oo 
(13) 
) 
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It can be demonstrated that the payoff function (14) or (15) represents the double output of a linear re- 
ceiver with the frequency response K(| f| — f,) if narrow-band noise with the spectral density s(| f | — f,) is 
fed to its input. 


Two particular cases of Eqs, (14) and (15) are given below: 
The case of broad-band noise 


@(A/)=S (df) K (Nas. (16) 


With an accuracy to the constant factor, the payoff function (16) reproduces the noise power spectral density, 


The case of narmww-band noise 


co 


® (Af) = K (Af) \ S(f)df. (17) 


Let us consider certain properties of the chosen payoff function. 


1, It is natural to assume that R(@ )H(c) in (12) tends to zero when the argument & tends to +a or —a@, 
Due to the fact that R(@ )H(c) is an even function, R(c )H(«) represents the Fourier transformation for the payoff 
function (12), Therefore, for Af — +00, the payoff function (12) tends to zero (perhaps while oscillating in the 
process), It is assumed here that the payoff function (6) also exhibits such an asymptotic behavior for large Af 
values, 


2. By definition, the payoff function is nonnegative and is an even function with respect to Af, i.e., 
(Af) = O(— Af) >0. 
3, It follows from the physical meaning of the payoff function that the first player (noise) must choose 


such carrier frequencies as to increase the values of the payoff function (5) or (12), and that the other player 


must strive to .cduce these values by a suitable choice of carrier frequencies, 


4, By differentiation with respect to f,; and f, in the f; = f, and f, = f, regions, it can be established 
that, if payoff function @(f,— f,) is nondecreasing with respect to f, in the f, = f, region, it is nonincreas- 
ing with respect to f, in the same region; in the f,; = f, region, the (f,;— f,) function is nonincreasing with 
respect to f, and nondecreasing with respect to f.. 


Therefore , by taking into account what has been said in paragraph 3, we arrive at the following method 
for manipulating the carrier frequencies, which is actually used in practice. 


While the carrier frequencies are being varied by the players, the noise must strive to reduce the differ- 
ence between its own carrier frequency and the useful signal-carrier frequency and, conversely, this difference 
must be increased for a better useful signal separation, ' 


Game Solution Methods 


The relationship between the payoff function analytical properties and the players’ optimal strategies is 
given in (7, 8]. 


Thus, if the payoff function consists of a “bell-shaped® function, which is, generally speaking, character- 
ized as positive, analytical, and monotonically decreasing to zero if its argument increases at a rate not lower 
than that of the exponent, the players" strategies are given by step functions [7]. 


The Gauss curve 


— fe) = exp [— (A — fa)? 
can serve as an example of such a payoff function. 
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Optimal strategies have also been found for payoff functions of the Green function type (8). Such func- 
tions are characterized by the fact that they are nonanalytical (it is usually assumed that there is a jump of the 
first derivative on the main diagonal of the square on which the payoff function is given, i.e., for f,;=f,). In 
this case, the players’ optimal strategies are formed from the continuous component that is different from zero 
throughout the entire [W,, W,] range, and the jumps at the W, and the W, boundaries, 


We shall now consider the application of general game solution methods [1,2] to the payoff functions (6) 
and (12), without imposing on the payoff function the special restrictions which were necessary in [7,8]. 


1, The theory of games provides the proof [1] that the solutions of the equations 


WwW, 
— ft (Add =, (as) 
w, 
\ — fa) (fa) df = (19) 


with the additional normalizing conditions 
Ww, 


Ww, 
\ (A) df, = \ (fs) df, = 1 (20) 


determine the optimal mixed strategies and the game value, i.e., 


= Fol), = l/r), 
in the case where 


> 9, >0 (21) 
for all f, and f2 contained in [W;, 


Equations (18) and (19) are Wiener-Hopf equations with finite integration limits and with a kernel which 
is an even function of the argument difference, Therefore, if the payoff function for the problem under considera- 
tion is approximated by the sum 


w 
® (41 — = exp(— — fal), (22) 
k= 


where ¢;, are constants and Rey, > 0 (k= 1,..., N), then, as is known, the solution of the Wiener-Hopf equa- 
tion is obtained in closed form [6]. After solving Eqs, (18) and (19), we must check whether condition (21) is 
satisfied for each concrete payoff function (22). 


Assume that the Fourier transformation of expression (22) is given by 


N—1, 


[] 
k=1 
N-1 N 
= 0, = N (23) 


s=1 k= 
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The solution of Eq. (18) [since Eqs. (18) and (19) are identical, the solutions will heneaiier be vainen only 
for (18)} will then be written in the following manner (6): 


= Ab (7, — W,) + + 
+ > {Cyexp [— — W1)] + Dj exp [— 15 (Wa — + Sp (/1), 

j=l 


where 5 is the Dirac delta-function, and ¢ ,,(f,) is the solution of the nonhomogeneous equation 


G[— (4) =». (25) 
For the case under consideration, this equation assumes the following form: 
t (26) 


By using, for instance, the method proposed in [6], it can be shown that the constants in (24) satisfy the 
following relations: 


A=B=av, C; = Dj = ew Gani, N —‘41), (27) 
where a and cj are determined by the poles and roots of G(r*), 
Thus, the final equations for the optimal mixed strategies of the players will be given by 


(/1) 18 —W,) + 8(7, + 


N-1 hen 
+ >} [exp | — — + exp 15 (Wa — + (28) 
j=1 
ns) = + Wal + 
N-1 
+ [exp | — — + expl + 5 (29) 
j=1 
The game value, determined from the normalizing condition (20), is equal to 
1 
j=1 
For instance, 
® (41 — fs) = (31) 
the optimal mixed strategies are 
Eo (/1) = (t+ — — Wa, 
(32) 
The game value is equal to 
29 


Let us briefly discuss the physical meaning of the obtained game solution, Since the solutions of Eqs. (18) 
and (19) usually contain delta-functions, it is more convenient to represent the optimal strategies in the form of 
a distribution function instead of in the form of probability density, For instance, the distribution function cor- 
responding to the probability density (28) is given by 


F = — Wi) + (f1 + at — Wel, (34) 


where 1[x] is the unit function, and c, and c, are constants which, in particular for the payoff function (31), are 
equal to 


2 + — Wi) 


Since the y parameter in (35) determines the "width" of the payoff function (31), we conclude from (35) 


that, if the payoff function width decreases, the probability portion for the jumps (in other words, the probability 
that the boundary frequencies W, and W, will be used) decreases and that the “value” of the carrier frequency 
uniform variation increases. 


Equation (33), which determines the game value, indicates that the noise effectiveness decreases as the 
y (Wz — W,) product value increases, Therefore, in the case of broadband noise, when its.spectral density coin- 
cides with the payoff function with accuracy to the factor [see (16)], white noise emission is the most effective 
one for noise, 


2. Carrier frequency changes within [W,, W,] usually occur in jumps, The consideration of such a game 
approximation leads to the so-called finite or matrix games, for which the solution methods have been explored 
to the greatest extent [2]. 


Let us divide the entire [W,, W2] range into (W, — W,)/h subranges, Then 
h=W,+ih, fp=W,+jh (0<i, 
In this case, the payoff function can be represented by a matrix whose elements a; are given by 
If the payoff function is different from zero only for 


<e, his (Wy, Wal, 
the game matrix is given by 


= for |i—j| (36) 
O(|i—j|h)=0 for |i—j| >=. (37) 


3. By using the S, N. Bernshtein polynomials for approximating* the payoff function, the solution of the 
game under consideration can be reduced to the solution of the so-called “separable” games, 


The payoff function for the separable game can be represented in the following form (2): 


f2) = Daw (41) 85 (Ja), (38) 


* The following theorem has been proved in the theory of games [2]. If the series 4(f.f2), O2(f1, fe), - + 
+++. ®,(f3, f2) of payoff functions that are continuous over the square, and which have game values equal to 
respectively, is uniformly convergent to the continuous,over-the-square payoff function ( f 
with the game value v , the vy, V2, ...,Un series converges tov, 


| 
| 
| | 


where a4; are constants, and 14(f,) and s;(f,) are continuous functions. The general methods for solving separable 
games are known (2). 


S. N, Bernshtein proved [9] that any continuous function f(x) that is defined in the 0 = x = 1 interval can 
be approximated by polynomials of the following form: 
= (39) 
while lim /,(z) = /(z), and Ci is the number of combinations of n with respect to m. 


As an example, we shall approximate the payoff functions (16) and (17) by S. N. Bernshtein polynomials, 
assuming that the noise spectral density or, correspondingly, the receiver frequency response, is a rational- fraction 
function. In this case, the approximating polynomial is given by 


= Car — | = ©, — = 
=0, Av) = () crave (1 — (40) 


where 0S v,51,05 v,=1, 
®,, (Av?) = = = (n), (41) 
i=0 k=0 


where a; (n) = (tym (2 


The general payoff function (6) can be reduced to the separable payoff function in a similar manner, 
Let us consider an example, Assume that 


(42) 


where \? (a7 wi is, according to the assumption, a small-value parameter. 
The S, N. Bernshtein polynomial for the payoff function (42) will be written thus: 


®, [(Av)*] = Cr Av™ (1 — = 
m=o 1+ —- 
k k n—m (43) 
=) (— 1)" m* Av — Avy. 
k=0 m=0 
It can be shown [9] that 


m* Ce Av (1 — Av?)"—™ = Av™ n* + 


m=0 


Thus, in the limiting case, we have 


lim ®, (Av*) = Av™. (44) 


k=0 


| 

- 


If the polynomial (43) is a second-degree polynomial (n = 2), for small values we obtain 
®, [(Av)"] = ©(0) + Av*[20 (+) — 20 (0)] 
+ Av! {© (0) + © (1) —20(+)} = 1 — (45) 
The solution of the game with the payoff function (45) is given by 


Mo (Vs) 8 (v2 —0) + — 1) = 


(46) 
+8 Wi) +8 (fs — 
The game value is equal to 
(47) 


The exact solution of the game with the payoff function (42) [2] yields the same optimal mixed strategies 


as the solution obtained from the approximate payoff function (45); however, the game value is equal to 


i 
v= wang (48) 
Further, the exact solution (46) is valid for 


Thus, the “approximate® game value is given by the first two terms of the series obtained by expanding 


the expression for the “exact” game value. 
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Certain numerical methods for determining periodical motions of automatic control systems are 
proposed here, For the application of these methods, it is sufficient to know the frequency char- 
acteristics of the linear elements and the diagrams of characteristics for the nonlinear elements. 


One of the main problems arising in theoretical investigations of the dynamics of automatic control systems 
is the determination and investigation of periodical motions, As a result of the efforts of many investigators ,con- 
siderable progress has been made toward solving this problem for systems that are close to linear systems (the 
small-value parameter method [1-3], the van der Pol method [4-7]), systems that satisfy the autoresonance or 
filter hypothesis (different variants of the describing function approach [8-11]), sectionally linear systems, the 
study of which can be reduced to the conversion of a straight line into another straight line [12-15], and for the 
simplest periodical motions of relay systems [16-19], However, in many cases encountered in practice, none of 
these methods can be applied. Moreover, even in those cases where the widely used describing function approach 
yields good results, there is some doubt about their correctness. 


The aim of this paper is to present a numerical method for determining periodical motions of automatic 
control systems which is based on the method proposed in [20]. The use of this method does not presuppose the 
knowledge of the analytical equations of motion for the system under investigation, It is sufficient to know the 
frequency responses of the linear elements, and to have the diagrams of the characteristics of the nonlinear ele- 
ments, 


1. Basic Initial Relationships 


Consider an open-loop system consisting of a linear element with the frequency response K(iw) = K,(w) + 
+ iK,(w), and a nonlinear functional element (Fig. 1), The frequency respouse K(iw) and the nonlinear element 
characteristic (x) can be given analytically as well as graphically. 


The periodical motions will be approximately assigned either in 

z Rs) y K(p) 8, the form of finite segments of a Fourier series, or in the form of a series 
of values at equidistant instants of time. 

Fig. 1. 


If we take into account n first harmonics, then, according to [20], 
between the Fourier coefficients A,, Az, ..., An, Bo, « « « » Bp of the periodical input quantity 


i= >} sin kot + By cos ket (1) 
k=0 
and the Fourier coefficients of the output quantity periodical component 
z=)>) A, sin kot + B, cos kot (2) 
k=0 
the following relationships hold: 
n 
al + +1, | x 
j=0 s=0 


x [K, (ko) sin — K, (ka}eos hj], 
j=o Le=o 
x [Ks (ke) cos + Ky (he) sin (3) 


where fo, fis++.» fm-y, are the values of the external input quantity f(t), which is applied together with x 
to the nonlinear element input at the instants of time 


Similarly, if x9, are the values of the periodical input quantity x, and Zp, 
are the values of the output quantity z periodical component at the equidistant instants of time 


tos by = by tm = fo + 
then 


n 
£5 = (ee) [Ks (ho) cos — 7) + Kel he) sin +h 
a=0 k=0 


In this, the number n of harmonics to be taken into account and the number m of equidistant instants of 
time in the 21/w period are chosen arbitrarily in correspondence with the required accuracy, but under the con- 
dition that m = 2n + 1 must hold, The m values equal to 12, 24, and 48 should be preferred if the calculations 
are to be simplified. 


An arbitrary automatic control system with a single nonlinear element can be obtained from the open- 
loop system (Fig. 1) by feeding the output signal z to the input, i.e., by assuming that z = x, 


Thus, in correspondence with relations (3) and (4), the equations for Ay, and By or xj, which approximately 
determine the periodical motion with the period 21/w in which we are interested, will be written thus: 


Ap= 9/54 sin + B, cos 


j=0 s=0 


x (hw) sin — K (kw) cos sil, 
(5) 


x (kw) cos + Ky (kw) sin 
and, correspondingly, also in the following form: 


m—1 n 
Q (2) (6 — j) + Ky (kw) k(s —f) +h= 


Equations (5) or (6) make it possible to find in the accepted appeonimetion the pestodical motions of the 
nonlinear system under consideration, 


(6) 


| | 


Actually, in the case of forced oscillations, the 21/w period is known, and Eqs. (5) or (6) constitute a sys- 
tem of 2n + 1 or, correspondingly, m equations with respect to 2n + 1 unknowns Ay, ..., Ap, Bo, «++» Bp OF, 
respectively, m unknowns Xo, ...,Xm-4- In the case of an autonomous system, i.e., if no external disturbances 
are present (f; = 0), the period 2x/w is unknown, The missing equation can be obtained by using the arbitrari- 
ness of the initial moment t, of time measurement, and by choosing it in such a manner that, for instanee, A, =0 
or B; = 0, For symmetric periodical motions, x, = 0 can be taken, 


The system of equations (6) can be written in the following form: 


m—1 
Fie (ao) + fis (7) 
o=Q 
where the matrix of coefficients aj, is given by 
Oy, Ag, Am—gy Sm—1 
Ao, My, ++ Am—g 
Am—1» Amy ++» (8) 
Og, 


and is, therefore, determined by its first row 


where 


n n an 

= K, (kw) cos ks + Ky (kw) sin™ ks = 2 ReS\K(ikw)e — (9) 
k=0 k=0 

if the nonlinear element characteristic, the external disturbance, and the periodical motion to be found 


are symmetric, then, for the even value m = 2m’, Eq. (6) will be written in the form given by (7) with the fol- 
lowing matrix: 


Om’—1» Oy, +» Om —2 
— Am’—2, — %m’—1, (10) 
— %, — Gq, — % 


where the values of & , are expressed according to Eq. (9) with the substitution of m' for m. 
Equations (5) can be written in the following form: 


2 
Ay + iB, = K (iko) ™ (11) 
j=o 
where 
= >} A.sin sj + B, cos —-sj + fj =Im2(A,+iB je ™ + (12) 
s=0 


For symmetric periodical motion and the even m = 2m’ value, the even harmonics vanish from Eq. (11), 
which is then written as 


m’—1 


Ay + iB, = =, K (iko) 
j=0 


3, (12°) 
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2. Method of Trial Hypotheses 


If the nonlinear element characteristic is composed of straight-line segments, the solution of Eqs, (5) or 
(6) can be found in the following manner, A certain assumption (trial hypothesis) is made about the position of 
the points Xo, » »Xm- On the linear sections of the nonlinear element characteristic, after which the 
system of Eqs, (5) or (6) becomes linear with respect to the Ay, and By or x; variables, respectively, This system 
is used for determining the unknown quantities, and the validity of the initial assumption is verified, If the ini- 
tial assumption is false, the next hypothesis is set forth, etc. 


Example 1,° Determine the forced oscillations of a proportional integral control system under the ac- 
tion of an external disturbance given by f(t) = 0,5sin 2t, which is fed to the input of a servomotor with the non- 
linear characteristic shown in Fig, 2, for the following parameter values: 


Ta=5, T,=0.01, B=1, 0.01, 5=0,1, @=1, Ty=100 


and a number of values of the T, parameter (T;* is the angular coefficient of the middle portion of the nonlinear 
element characteristic), The control system structural diagram is shown in Fig. 3. In the case under consideration, 
P(T, P+ 8)(T, p+ + +9) 
The transfer locus K(iw) is shown in Fig. 4. We shall use the system (5), assuming that m = 12 and n = 5, 
By taking into account the symmetry of the periodical motion to be found, this system will be written as 


= 10.T40 (24) + (23) + 0,602 (24) + 0,260 (24) — 0,20(24) —0.290 


= 50.202 (a) + 0.740 (21) + (2s) + 0.60 (a) + 0.260 (x4) — 0.20 + 0.25, 

= 10.20 (0) + 0.290 (21) + 0.740 (a) + (20) + 0.602 (x4) -+ 0.260 + 0.42, (14) 
[— 0.260 (25) + 0.20 (a3) + 0.290 + 0.740 (as) + + 0.60 + 0.5, 

[— 0.682 — 0.260 (21) + 0.20 (23) + 0.2982 (zs) + 0.740 (x4) + + 0.42, 


= [— (24) — 0.682 (21) — 0.280 (x4) + 0.20 (x3) + 0.200 (x4) + 0.740 40.25, 

The remaining values of,x; will be found from the equation xj = —xj.. (j = 6, 7,8, 9, 10, 11, 12), 
Let us consider the following hypotheses: 
a. for T, = 0.45, 

Q (x0) = Q (x1) = Q (x3) = Q (zs) = 1, Q(x) = = Fm, 
b. for T, = 0.33, 

c. for T, = 0.1, 

Q (x) = (21) = (x2) = (zs) = Qs) =1, 

After substituting the (x;) values in the system of equations (14), we find that, for the hypothesis a, 


= 0.95; 2, = 1.11; 2g = 1.03; zy 0.58; 
—O0,15; 2 = —0.2; 


* This and the following examples were composed by A, Al'tman and M, Yudkovich, 
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for the hypothesis b, 
= 0.99; 2,=—1.16; zg=—1.0; 
m=—O01; =—0.5; 
and for the hypothesis c, 


= 0.89; 2, = 1,28; = 1.27; = 0.96; 
a=0.21; = —0.45. 


In each of the a, b, and c cases, the hypothesis used proved to be correct (the rejected variants are not 
given here), 
Similar calculations for m = 24 and n = 5 yielded the following results: 
a. for T, = 0.45, 
= 0.92; = 1,12; = 1,2; a =0.93; 25 = 0.76, 
0.5; 27 =0.2; = —0.24; = —0.4; = — 0.66; 
b. for T,; = 0.33, 


= 1.0; %=1.07; <2 =1.10; 2=1.05; 25 = 0.70; 
= 0,4; 2, = 0.08; = —0.39; x19 = — 0.56; = — 0.72; 


Cc. for T; =0,1, 


= 0,94; 2, =:1,13; =—1,25; = 1,27; = 1,2; 25 = 1,06; 
= 0,86; 27=0,5; 2e=—0,11; z——0,2; 29 = —0,46; = — 0,67. 
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Example 2, Consider the same proportional integral control system in the absence of external disturb- 
ances. The equilibrium state of the system under consideration is stable for T, > 0,59, and it is unstable for 


T,.< 0.59. Let m = 12 and n = 5; then, if we take into account the symmetry of the periodical motion to be de- 
termined, the matrix of the system of equations (5) will have the form given by (10), where 


1 
= [Ki (@) + Ki (30) + Ki (5e)}, 
a1 [0.85K, (0) — 0.85K; (50) + 0.5K, (0) (80) + 0.5K, (5a)}, 
[0.5K (@) — Ky (3@) + 0.5K 1(5@) + 0.85 Ke (@) — 0.85K 


(15) 
a3 = (@)—K, (3@) + Ky (5@)}, 


1 
= [— 0.5K (w) + Ki (30) — 0.5K; (50) + 0.85Ks (@) — 0.85K (5e)], 


1 
= [— 0.85K, (@) + 0.85Kj (5@) + (@) + Ka (30) + 0.5K (Se)}. 


In addition to these equations, we shall assume that xp = 0, Let T,; = 0.45 and T, = 0,33, For the first case, 
we shall assume that 


Q(a)= Q = Q (zs) = Q(x) = 1, Q (es) = 


In correspondence with this hypothesis, the system of equations for x4, Xp, Xs» Xq» Xg, and w will assume the 
following form: 


= +. ag + Og + + = + + + + 
= — + Oo + + + 2.20925, = — — Og + Ay + + 2.20475, 
= —2.2agry — — Os + + 2.225, = — — ay — — + 


where the Xs, and coefficients are given by Eq. (15), 


If we express the unknowns x, and xg from the first and the sixth equation in terms of w, and substitute 
them in the second equation, by using the graphical method, we find that w = 2,6, For this value of w, we find 
Xg, Xs, and x, from the remaining equations, Thus, 


%=0; =0.34; zy =—0.68; =—0.61; 0.4 


and the hypothesis used has been proved correct, 
For T, = 0.33, and assuming that 


Q (21) = Q (x2) = Q (zs) = Q (xs) = Q (ws) = 1, 
we obtain the following system: 


O = + Og + Os, = + + + + 


From the first equation, we find that w = 2,6, and, from the remaining equations, we find that 
a =0,38; 23=0.55; 23=—0.73; 2%=—0.63; 2 0,35. 
Similar calculations for m = 24 and n = 5 for the case where T, = 0.45 and T, = 0,33 yield, respectively, 


@=2.6; 2 = 0.15; 0.36; 23 = 0,48; = 0,57; 25 = 0.62; 
%e=0.7; 2 =0.69; 2 =0.5; 0.31; = 0.2 


and w=2.6; 2 =0.2; = 0.38; 
2y=0.47; =—0.59; = 0.68; 

%=0.73; =0.72; = 0.66; 

=0.54; 29 =0.34; 2, = 0.17. 
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For T, = 0, the nonlinear element becomes a relay element, and the periodical regime can be found ex- 
actly [16-19], For the sake of comparison, Fig. 5 shows the diagrams of the exact and the approximate solutions, 
The approximate solution (curve with white points) was found for m = 24 and n = 5, 


3. Method of Consecutive Approximations (Iteration) 


The form of equations (5) and (6) is suitable for the direct application of the iteration method, In the case 
of a nonautonomous system, we take any values of Aj and Bj or xj for the initial zero approximation and, by sub- 
stituting chem in the right-hand sides of equations (5) or (6), respectively, we find the first approximation; in the 
same way we find the second approximation by using the first approximation, etc. If the thus obtained consecu- 
tive values converge, they converge toward the solution, For an autonomous system, the application of the itera- 
tion method is somewhat more complicated; it can be done in the following manner, We write system (6), to- 
gether with the additional requirement given by 


= (@) + F (2g) 2) +++ Lm—1; @) = 0. (16) 


For the assigned zero approximation for x», X;, «. +» Xm-4 and w, we find the values of x; from the first 
group of equations by using the iteration method in the manner described above; then, by substituting the found 
values of x; in the last equation, we find the new value of w, From the first group of equations we find the next 
approximation for the x; variables, etc, 


For the Ay and By, variables, the iteration method can be applied in a similar manner. In the absence of 
convergence in the iteration of all Aj, and By, variables, it is possible to secure convergence in a number of cases 
by iteration only with respect to the Fourier coefficients of the higher-order harmonics, 


Example 3, We shall determine the self-oscillations for the same control system as in the previous 
example with the difference , however, that the nonlinear element has the characteristic shown in Fig, 6, For 
m = 12,n = 5, and x, = 0, we arrive at the following system: 


O = (0) + (ay) + (x2) + (xs) + (x4) + (x5), 
= — + ao (21) (2g) + (xs) + (x4) + (x5), 


= — a,Q (0) — (2) — (72) — (as) — (x4) + (25), 


where values of %» to &; are determined by Eq. (15), For the initial approximation, we shall take x, = x, =...= 
= Xs = 0.5. Then, from the first equation, we find w = 2.6 and, from the remaining equations, 


a, = 0.33; 2 = 0.57; x3 = 0.72; x = 0.67; 25 = 0,38. 


For the new values of x, X2, . . « » Xs, We find w = 2,62 from the first equation and, from the remaining 
five equations, we find x, = 0.31, x, = 0.51, xs = 0.67, x, = 0.61, and xg = 0,35, 


The third approximation yields x, = 0.31, x, = 0.54, xs = 0.66, x, = 0,60, and x, = 0,35, 


The obtained solution can be rendered more accurate by using m = 24 and n = 5, and by taking the just- 
found values for the initial approximation. The intermediate values are found by linear interpolation, After two 
iterations, we obtain: x, = 0.13, x, = 0.28, xs = 0.41, x, = 0.5, x5 = 0,58, xg = 0.61, x7 = 0.61, xg = 0.54, xX = 
= 0.46, Xgq = 0.33, and xq, = 0,20, 
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Example 4, Assume that the nonlinear element has the asymmetric characteristic shown in Fig. 7. In 
system (5), we shall put m = 24 and n = 5, and we shall additfonally assume that B, = 0, Since K(0) = ow, By = 0, 
By taking Ay = 0,34, Ap = Ay =... = As = 0, By =. ~~ = Bs = 0 for the initial (zero) approximation, and by using 
Eqs. (12") and the graph in Fig, 7, we find, in the order of increasing values of the number j , the values of O(x;): 
0, 0.3, 0.53, 0,72, 0.83, 0.88, 0.9, 0.88, 0.83, 0.72, 0,53, 0.3, 0.0, -0,18, -0.24, -0,28, -0.31, -0.33, -0.34, 
—0,33,—-0.31, -0,28, -0.24, and -0,18, 


By substituting these values of (xj) in the equation B, = 0, we find that w = 2,6, After substituting the 
found value of w in Eq. (5), we find the following approximation for the Ay, and By Fourier coefficients : A, = 
= 0.39, Ap = 0,02, Ay = Ag = Ag = 0; By = B, = 0, By = 0.01, By = Bg = Bg = 0, 


The third approximation is a repetition of the foregoing with the accepted accuracy; in this, w = 2,62, 
Ay = 0.42, Ap = 0,02, Ay = Ay = Ag = 0; By = By = 0, By = 0,01, By = 0.01, By = Bg = 0, 


4. Method of Extension with Respect to a Parameter 


The extension of the solution with respect to a parameter is a y2ry flexible and efficient method for de- 
termining the solutions of equations (5) or (6), In this, che parameter with respect to which the solution is ex- 
tended can be either one of the parameter: of the system under consideratica or a parameter artificially intro- 
duced for extending the solutions from equations with known solutions to equations whose solutions are to be de- 
termined, 


The combination of the above-described methods with the extension with respect to a parameter is con- 


venient, since the required hypothesis can always be readily found, or a good initial approximation is already 
available, 


A more independent application of the method of extension with respect to a parameter can also be used, 
In the general case, if we assume that K(iw) and (x) are functions of a parameter !, and if we differentiate 
the system of equations (5) or (6) with respect to this parameter, we obtain the following system of differential 
equations: 


is 2% 
ie 2% (j—2) 
+ Q.( ste) Re |isKy (isw) + (iso) Je" ™ } (7) 
or, correspondingly, 
m—1 an 
where 


Assume that, for By using this solution 
for the initial conditions for |! = [l9, and by solving the system of differential equations (17) or (18) by one of 
the approximate methods, we find the necessary solution for a certain interval of the parameter }! values, 


It should be noted that, in the case where the nonlinear element has a sectionally linear characteristic for 
all variables except w, the systems of equations (17) and (18) represent systems of linear differential equations 
with constant coefficients for each extension interval where no changes occur in the arrangement of Xo, X;, + «+ 
Xr Values on the linear portions of the nonlinear characteristic, 


In conclusion, it should be emphasized that the application of the method of extension with respect to a 
Parameter is connected with studies of possible bifurcations of periodic motions [21]. 
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DETERMINATION OF AIRCRAFT CONTROL EQUATIONS FOR THE 
OPTIMUM PATH UNDER VARIABLE-WIND FLIGHT CONDITIONS 


V. D. Matytsin. 

Translated from Avtomatika i Telemekhanika, Vol, 22, No. 1, 
pp. 57-66, January, 1961 

Original article submitted April 20, 1960 


The present article is concerned with the extremal problem in determining the conditions neces- 
sary for the quickest aircraft flight from one point in space to another, under conditions where the 
wind changes with respect to coordinates and in time, The minimum flight time is taken as the 
criterion for the optimum path, Variation calculus is used as the mathematical tool (the relative 
extremum variation problem), 


If the effect of wind along the flight path is given detailed consideration, the advantage of the shortest- 
distance flight becomes doubtful, It is sometimes more advantageous somewhat to lengthen the route if, by tak- 
ing advantage of the wind, the over-all flight duration can be reduced, 


The determination of the optimum flight paths is of practical interest from the point of minimizing the 
flight time. The knowledge of such flight paths makes it possible to obtain more complete information on the 
maximum aircraft range, and to perform the planned flight with minimum fuel consumption [1]. 


At the present time, cruising flights are conducted at altitudes where the wind velocity can attain con- 
siderable magnitudes, Consequently, the equations of motion must reflect, to the fullest possible extent, the ef- 
fect of wind on the aircraft dynamics, 


In connection with the fact that flights take place practically at the same altitude, we shall consider the 
aircraft motion problem in the course plane, The coordinate 
systems and the basic notation (see figure), which are used in 
correspondence with we accepted standard [2], are given below, 


The following coordinate systems are used: XOZ is the 
terrestrial coordinate system, X,OZ, is the bound coordinate sys- 
tem, and X,OZ, is the velocity coordinate system, 


In the X,0Z, coordinate system, the OX, axis has the di- 
rection of the flight velocity vector, 


In composing the equations of motion, the following nota- 
tion was used: V, is the path (absolute) aircraft velocity, V, is 
the aircraft air (relative) velocity, W is the wind velocity, P is 
the aircraft propelling force, Q, is the front drag force, Z, is the 
lateral force, ¥ is the aircraft yaw angle, ¢, is the direction of 
the air velocity vector, ¢ is the aircraft course angle (direction 
of the path velocity vector), 8, is the relative (aerodynamic) 
slip angle, which is formed by the positive directions of the OX, 
and OX, axes, 8, is the absolute slip angle, which is formed by the positive directions of the OX, and OX, axes, 
and 6 y is the additional slip angle due to the wind, which is formed by the positive directions of the OX, and 
OX; axes, 


The aerodynamic forces and moments depend on changes in the air (relative) velocity; they are given by 
the following equations: 
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ta 
W, 

\ | 

= 


Q, = V3 
Z, = C, = (1) 


l 
M, = m,* iv? = y? [ms Br + mi + | 
Under actual flight conditions, the wind exerts the most frequent and the longest effects on the aircraft. 
The wind can be defined in a terrestrial coordinate system by its projections® : 


W =({W, (a, z, t), 2, 


In the aircraft motion equations, we shall take into account the fact that the wind changes the aerodynamic 


forces acting on the aircraft by changing the aircraft velocity with respect to air, without, however, affecting the 
inertial terms of the equation of motion. 


The inertial terms depend on the terrestrial (absolute) velocity and acceleration, and the aerodynamic 
forces depend on the air (relative) velocity. 


S. B, Puzrin was the first to derive the aircraft motion equations by taking into account the wind effect [3]. 
We shall use his method, 


Let us project all forces acting on the aircraft on the OX, and OZ, axes, Projection on the OX, axis yields 


m = P cos — Q, £08 By-+ Zr sin By, (2) 
and projection on the OZ, axis gives us 
mV, = — Psin®, + Q,sinB,+ Z, cos By (3) 


where Rg is the instantaneous curvature radius of the absolute path, 
The equation of moments with respect to the y axis has the following form: 


ly st = My. (4) 


The system of Eqs, (1)-(3) is supplemented by the kinematic equations 


= W, + V,cos@, = V,cos a, z= W,—V,sing, = —V,sing,, 
5 


Or = Bata, Bw = Ba — Br = Gr — Pa- 
By using system (5), we find the expression for the path velocity: 
=Vwi4+w?+ V? + 2V,(W;.cos 9, — W, sin 9r); (6) 


Considering W,/V, and W,/V, as small quantities of the first order, we obtain, by transforming expression (6), 


> v2 


r 


sv, (14 


r 


* Flight in the horizontal plane is considered. 


m (6), 


Thus, 


V.=V,+ W, — W, sing. (7) 


System (5) also makes it possible to determine the additional slip angle By, which is due to the wind: 


V, cos (@r — Bw) = Wx + V,cos @, = V, cos @, cos By+ V, sin @, sin By . 
(8 


— V, sin (9; — Bw) = W, —V,sing, = —V, sing, cos By+ V, cos @, sin B,, 


By multiplying the first equation in system (8) by sing,, and the second equation by cos¢,, we obtain, 
after addition: 


V,sinBy= W,sing, + W, cos@, 


W,,sin 9, + W, cos 9, (9) 
V . 


sin 
a 
By taking into account (1)-(7), aud the smallness of angles By and 8;, we have 
V,=V,+ W, cos — W,sing,, 
dp, Bw r 
mV, ( = — + Bl + QrBw-+ Zr Br. 
v2 
= + my 8) + SPL 
By replacing the equation of motion about the center of gravity by the balancing relation, we finally obtain 
V,=V,+ W,.cos 9, —W,sing,, 


4B 
mV — = (— P+ + (—P + 
Ww 


v a 
z= W,+V,cos@,, t= W,—V,sing,, 


Assume that the aircraft is at the point A(x», Z) at the initial moment of time, and that the flight destina~- 
tion is B(x,, z,). The aircraft flies under conditions where the wind velocity depends on coordinates and the time, 


We shall try to determine the rudder deflection equation, which will make it possible to complete the 
flight along a path passing through points A(x», Zo) and B(x, 2,) in the shortest possible time, The existence of 
such a path follows from the physical meaning of the problem, 


The equations of system (10) are differential linkage equations, The analysis obviously cannot be com- 
pleted for the general case, due to the complexity of the linkage equations, Therefore, we shall solve this pro- 
blem in two stages, Considering ..ds problem as a purely kinematic problem, we shall determine the equation 
of changes in the angle ¢, of the air velocity vector in order to obtain the optimum path, Then, by using the 
obtained dependence and the dynamic equations, and by taking into account the wind, we shall determine the 
optimum rudder deflection equation, 


Considering that the air velocity vector value does not change considerably in flight, we shail consider V, 
as a Constant quantity, For these conditions, the problem was solved by I, M, Gyunter [4], who reduced it to the 
so-called Meyer variation problem, We shall now give a simpler and shorter derivation of the equation of changes 
in the air velocity vector angle for the optimum path, 


We shall formally consider the flight time as a quantity depending on a certain parameter € , t = t(€ ); then 
x = x(€ ), z = z(€ ); at the initial instant of time, 


(10) 
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= (Eo), =2(Eo). =2(Eo), 
and at the final instant of time, 


The time differential will be given by dt = “dg , and the flight time will be determined by the integral 


dt (E) dk. (12) 


The kinematic linkage equations will assume the following form: 


=; (Ws + V,cos Gr), de —V,sin r)- (13) 


Following M. A. Lavrent'ev and L, A. Lyusternik [5], we shall denote the dt/d€ function by F, i.e., 


Equation (12) will assume the following form: 


Ey 
FUE, 1, 2, Les Gre) (14) 


and the linkage equations (13) will become 


= (Wz + V,cos Fy = (W,—V,sin F. (15) 


The flight time (14) can be considered as a functional of F. The F function itself will be determined from 
the first equation of system (15), namely 


16 
P= 
The flight path in which we are interested must satisfy the linkage equation which will be obtained by 
eliminating the F function from (15): 


wWw,—V,sing 
x, 2, Te, ze) = — =0. (17) 


The class of the acceptable paths to be determined will consist of curves satisfying the differential linkage 
equation (17) and the boundary conditions (11), We shall now find the condition for the extremum of functional 
(14). For the solution of this problem, we shall use the following theorem, 


The x, z, and ¢, functions that secure the extremum of the functional 


Fe, ©, 2, Pry Ley Pre) 


in the presence of the linkage equation 


2, Pr, Ze, Ze, Pre) =0 
satisfy for a suitable choice of the \(€ ) factor the Euler equations composed for the functional 


| 


The A(E), x, z, and ¢, functions are determined from the following Euler equations; 


0H d 0H d OH 
and from the linkage equation*® 
W,—V,sin 


Wz +V, 08g, 
In our case, the system of equations (18) assumes the following form: 


On, 


It is essential that A be a function of € only. 


By using expression (16) for F, we shall determine the 0F/ 0x, OF/0z, OF/ d¢,, OF /OK, , OF /O%, , and 
OF /89,, derivatives: 


aw, 29, 09, 
— V,sin ¢, =— V, sin 9, 
dz (W, +V, cos W,+V, cos 9, 
ow r ow, 39, 
aF —V, sin 9, F V, sin 9, 
39, * +V,cosq,” We+V,cos@, 
OF i OF OF 


Similarly, by using the linkage equation (17), we obtain 


= V,sin — F (3! —V, cong, 


W,—V,sin 9, OW, 39, ow, 09, 

—V sin — F (“Gt — Veco (21) 
W,—V,sin 9, 
— FV,sin 9, + FV,cos 


By taking into account (7), the third equation of (21) can be reduced to the form 


* See [5], Ch, V, paragraph 22; also [6], Ch, IV, paragraph 2, 
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Since all the necessary mathematical operations for determining the optimum law of ¢, variation have 
been performed, we can put t = € in systems (20) and (21), Then, = = In this case, it 
follows from (10) and (16) that F = 1, 


For F = 1 and by taking into account (15), (20), and (21), system (19) will have the following form: 


ow 
sing, W,—V,sing, (Ws y 
W, cos 9, 9, ( Ox Ox )- 
—1(54—V,cos 9, +V, cos @, —h | = 

ow ag 

sin @, W,—V,sing, (Ss 

ow, 99, dh 
GE) — =0, 


V,.sin @, 
+V, cos @, W. +V, cos 9, 


z=W, + V, cos @,, z= W,—V,sing,. 
We shall determine \ from the third equation of system (22): 


sin 
(23) 


The first equations of system (22) can be written as 
ow. Vv 09, [ 1 W,—V,sin 
ow, 09, d 1 W,—V,sin 9, 
— | — hw | = 0, 
ow, 09, W,—V,sin@, 


—V,cosq, =0. (24) 


By using (7) and (23), we find that 
4 W,—V,sing, W,, cos 9, —V, sin’ 9, 
+V,, cos 9, W,, +V,, cos 9, V, (W, +V,cos 9,) 
cos* 9, + W,, cos 9, COS (25) 
~ 
System (24) assumes the final form: 


ow 
— 9, +V, sin cosy, + sin —V, sin @, cos @, —— 
V 


aw ap, ow 
V 
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ow ow, 
— sin 9,9, + cos 9, — sing, 


(26) 


Ow, aw, 
£08 @,9, + sin — COs @, 0 


By multiplying the first of the equations (26) by sing,, the second by cos¢,, and by subtracting the first 
from the second, we obtain 


ow, ow aw aw 
+( — ) sin + sin® cos" = 0. (27) 


It should be noted that the difference between the signs in Eq, (27) and those in the results obtained in [4] 
) can be explained by the different choice of the z axis positive direction and, correspondingly, the positive direc- 
tion of the W, wind component. The coordinate axes used in the present paper are those recommended by the 
standard [2]. 


The above derivation is shorter, which is due to the simpler form of the functional under investigation , 


F= 
W,+V,cos@, 


in comparison with 
[V2 -+ W2 + F? + 2(Wate + Wy F — (22 + yt) = 0, 
that was used in [4]. 


Thus, over the optimum path, the air velocity vector changes its direction according to Eq, (27), 
Finally, for determining the rudder deflection equation, we have the system of Eqs, (10) and (27): 


) 


W, cos @, + W, sin 9, 


| 


Va = V,+ W,cos 9, — W, sing,, Bw= 


mf" 
ma — =(—P + + (—P+Q)By, 
v 
+ ) sing, cos + — cost = 0, (28) 


z= W,,+V,cosq@,, z= W,—V,sing,. 
We shall find By from the second of the equations (28) 
gs W,,V, sing, + 9, +W, Vacos 9, — W,V, sin 9, 9, 
y2 


_ sing, + W,V, 
Va 


or 
ow ow 
aw, aw, 
or 
= 
| 
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or, by expanding V,, we obtain 


W,csg,—W,sing, W,.V, sing, + W,V, cos @, 


W,, sin 9, (W,, cos @, — W, sing, 9, — W, sing, — W,, cos 9, 9,) 
W,, cos 9, (W,, cos — W,, sin — W, sin p, — W, cos 
By using Eq. (7), 
W,,cos 9, —W, sing, =V, —V,; 
after regrouping, we obtain 
a 
4 


a 
W,(V, cos 9, + W, sin* p, + W, sin @, cos @,) 


or, by substituting Va= Vr + WxCos¢r— Wzsin¢,, 


V,V, —V,) + (W,, sin 9, + W, cos 

Bw v3 + 

(V, sin 9, + W,, cos @, sin 9, — W, sin® — W,, sin @, cos — W,, cos* 


a 
W, V, cos @, + W,, cos — W, sin cos + W,, sin? p, + W, sin @, cos @,) 
= > 


Neglecting the higher-order terms, W3/V?, W,W,/ Vi, after regrouping, we have 


+ 


Va—v,. — WwW, (29) 
Bw= Gr + y2 


Bow r= Or 


By using the second and the third equation of system (28), as well as (30), we obtain. 


(30) 


| W,—2W, —P+2," (— P + Q,)(W,, sin + W, cos @,) 
Br mv? 


WA—W,2 (—P+Z5")B, (—P+Q,)(W, sing, + W, 00s 


By taking into account the fourth and the fifth equation of system (28), the equation for determining the 
optimum rudder deflection law will have the following form: 


_ 
Q,) (W, sing, + W, 
mV,V, 


By expanding W, and W, , after regrouping, we finally obtain the rudder deflection equation: 


—P+Q, W,sing, + W, cos 9, 
my Ma 
mir m 


+ V,cos — (Ws — V,sin,) + 


aw, aw, 
( Oz 


24. W24V,W,cos*@, + V,W,sin* 9, — 3V,W,sing,) + 


W,+ W,V,cos*® 9, —V,W, sin® -) — 


aw 
+ + Wi—V,W, —V, W cos? g, + 3V, 


where ¢, and V, satisfy the equations 


ow ow 
V,=V,+ W,,cos9, — W,sing,. 
The optimum path coordinates are determined by the equations 


z=W,+V,cosg,, z=W,—V,sing,. 
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AUTOMATIC OPTIMIZATION OF SPATIAL DISTRIBUTION.® I 


L. N. Fitsner 

Translated from Avtématika i Telemekhanika, Vol, 22, No. 1, 
pp. 67-76, January, 1961 

Original article submitted July 22, 1960 


The principles of automating the spatial distribution of some physical substance are examined. 
Certain theoretical questions concerned with the solution of the stated problem by means of auto- 
matic search are discussed, 


1. Class of Solvable Problems 


It becomes necessary in solving some technical problems to distribute some physical substance @ in three- 
dimensional space in such a manner that the values of @ will approximate some requirement at n given points in 
the space in the best manner, 


For example, let it be necessary to distribute the material allowed for machining parts of complex shape 
from blanks. To do this, the blank is moved about in space so that the individual points on its surface approxi- 
mate the corresponding points in a space given beforehand in the best manner, Here the distance between the 
points on the surface of the blank and the corresponding ones in the space is optimized, The problem of distri- 
buting flows in space is quite close to the given problem in the method of its solution, 


Let us assume that there are several sound sources located at m points in the space, and that it is necessary 
to find the power of each of the sources in such a way that the sound energy at n other points approximates cer- 
tain required values in the best manner. 


Analogous problems may arise in distributing the flows of other physical things, such as light, radio waves, etc, 


The principles of automation discussed below can also be applied in automatic parts assembiage, in the 
pattern layout of materials for stamping, etc, 


We notice certain peculiarities in the class of technical problems being examined, The values of © at in- 
dividual points in the space are interrelated and, as a rule, cannot be set independently, The number of points m 
at which the values of @ are given, and the number of mechanisms n which act on @, are not equal in the general 
case, Ordinarily, it is not possible to get values of 6 which precisely agree with requirements if'a limited num- 
ber of points are available; as a result, the problem solution is only an optimum approximation from one point of view. 


It is assumed that the required distribution law does not change with time, 
and that the problem is considered solved when the best approximation of @ to 
the given distribution is attained, Deviations of @ from the required values can 
be measured at all m points of the given o. 


All the results of the study presented below are related only to one pro- 
blem, important from a technicoeconomic view , and concerned with the distri- 
Fig. 1. bution of the material allowed for machining parts of complex shape from blanks, 
However, a large part of the results obtained are applicable also to other problems of 
this class, 


2. Distribution of Cut Material 


In machining the blanks of parts on a machine tool, it is important to position them in the best manner, so 
that the layer of material removed will be approximately identical at all places on the surface being machined. 


* This article was reported at the IAT of the Academy of Sciences, USSR, in the Session on the Theory of Opti- 
mum Self-adjusting and Mutually Interconnected Control Systems, Automatic Control Theory Seminar, May 11, 
1960, 
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This will result in the smallest amount of material to be cut away, thus permitting a reduction in the machining 
time and amount of waste, If the blank has a complex shape, it is extremely difficult, as a rule, to chuck it in 
the tool; many auxiliary attachments are necessary, and even with their aid painstaking labor is required, There- 
fore , it is important to automate this process, 


The outline of a blank (solid line) and outline of the part (broken line) are shown in Fig. 1 (for simplicity, 
a flat part will be considered), Let the blank be moved without deformation in three degrees of freedom (pro- 
gressive movement along x, and x2, and rotation x, around an axis perpendicular to the plane of the figure), The 
problem is to vary xX, X_, and xs so that the best, from the viewpoint of some criterion, mutual disposition of the 
broken and solid outlines will be achieved, 


Let us assume that removing a uniform layer of cut material results in an ideal part, We will call such a 
blank outline an isoform outline of the ideal part, In this case the blank can be set so that the normals AB and 
A,B,, erected at any two points on the ideal part outline, will be equal to one another, Such a mutual disposition 
will be ideal, However, in general, the blank outline is not an isoform of the ideal part, and the problem can 
only be one of finding the best approximation to the ideal case, from the standpoint of some criterion or other, 


Let us consider three criteria of approximation having practical interest. 


Let us designate the distance (AB) between the outlines by €;, and assume that measurements are made at 
‘m points along the part outline. We will assume that €; is always positive, Then the first criterion will be 


max %) (i =1,2,...,m). (1) 
The optimum positioning problem consists of choosing x4, Xz, and x3 so that Q, will be minimized: 


minQ, = Ei (2, 2, (i =1,2,..., m). (2) 
Xt, 


In the second criterion, the minima of €; are maximized, This kind of blank positioning eliminates the 
dangerous condition of having negative cut material: 


Qs = min = min (%, (i = 4... , m). (3) 
It is necessary to find x4, X,, and x, which maximize Q,: 
max Q,= max {min§;(x, 2%, %)}  (i=1,...,m). (4) 
(X1, Xa» Xs) (X14, Xs) 
Finally, the third criterion may be the dispersion of £;, i.e., the difference between the maximum and 
minimum: 


Qs = Qe = max bi (x, Te, — min Zs) (5) 


(i= i,...,m), (fj =i,...,m). 


If this criterion is used, the optimum positioning problem consists of minimizing it: 


min Q;= min {max (x, — min 2, %)} (6) 
Xa, Xs) Xa Xs) 


(é=i,...,m), t,..., me). 
If the part blank is isoform, all three criteria result in a uniform cut, All the considerations presented re- 
late to the case where three measurements are made, 


3. Automation Principles 


The machining process can be automated with three types of automatic control systems: systems without 
feedback, automatic control systems, and automatic search systems, We will deal with each of these systems, 


Figure 2 shows two isoform parts: one by a solid line (A), and the other by a broken line (B), The distance 
between the sides of pentagons A and B can be measured by the indicators  ;~€5, and part A can be moved along 
the axes x, and x, by systems without feedback, The displacement along the axes mentioned can be found from 
the following formulas: 
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Fig. 3, 


Az, = —4, Az, = — 6. (7) 


It is assumed here that the size of cut material 5 is known beforehand, This limitation is not necessary if 
two automatic control systems are used for positioning, the first of which moves the part along x, into a position 
where 

=0, (8) 
and the second determines the position along axis x,: 


(9) 


Let us examine a more general case of optimizing the positioning of isoform parts, If the position devia- 
tions from the desired optimum are small, and if the functions 


Ei Xa, » (i=1,...,m) (10) 
have continuous partial first derivatives, the following set of equations is justified: 


+ 


The coefficients of 0€4/0xj i=1,...,m;j=1,...,n) are assumed to be constant and can be deter- 
mined beforehand experimentally for the given set of equations, The set of algebraic equations (11) can be solved 
by one of the known methods on a computer [1], and the values of x; found can be used to govern the actuating 
mechanisms of the system by means of n systems without feedback, or n completely independent automatic con- 
trol systems, Data on the cut material size is not needed if this method is used; isoformity of the parts and 
stability of the coefficients of 0€ ; /0x; in time are necessary. 


Automatic search is the best general method for solving the given problem, and is useful for the case of 
optimizing the machining of nonisoform parts, A block diagram of a system for optimizing this process by the 
automatic search method is shown in Fig. 3. Platform P with part D fastened on it is moved in n degrees of free- 
dom by actuating mechanisms 10,, 10,,...,IOn. The deviation of individual points on the blank from a posi- 
tion given by a template is measured by indicators I,, I,,...,Im. The outputs of the latter act on the input of 
the optimum positioning criterion block (OPCB), which generates a signal Q in correspondence with one of the 
criteria (1), (3), or (5), The automatic optimizer AO excites the actuating mechanisms (10,, I0,,... .10Op), 
and the minimums of Q, or Qs, or the maximum of Q, is found by means of automatic search, 
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4, Theoretical Questions Related to the Automation 


of the Spatial Positioning 


A number of theoretical questions arise in solving problems related to the optimization of spatial position- 
ing by the automatic search method, The main ones are the following: a) choosing the method of automatic 
search; b) evaluating the effect of the system inertial elements, and also the static characteristic of the search 
object on the stability and accuracy of the search process; c) evaluating the effect of backlash and stick-slip 
friction; d) evaluating the effect of random disturbances; e) evaluating the effect of the insensitive zones of 
the optimizer; f) finding the optimum form of trial probes; g) synthesizing a coordinate converter which de- 
creases the search time; h) finding the optimum scan amplitude and period; and j) finding a search law which 
is optimum in use of time, 


The greater part of these problems can be studied analytically, However, it is not desirable to use analytic 


methods for investigating some of these problems, since it may be necessary to develop special apparatus to synth- 
esize automatic search systems automatically; these devices are called automatic synthesizers, 


The theoretical studies of optimum positioning of automatic search systems which have been enumerated, 
done both analytically and with the aid of automatic synthesizers, and also a description of the operating prin- 
ciple of the synthesizer itself, will be the subject of succeeding articles, 


5. Concerning the Choice of the Automatic Search Method 


It is difficult to base the choise of search method for the optimum positioning of blanks during machining 
on a strictly mathe matical analysis, and we will therefore limit ourselves only to general engineering considera- 


The accuracy and time of search, and the com- 
plexity of the actual control apparatus required, are 
ordinarily the basic choice criteria, A number of 
the specific peculiarities of the system must be con- 
sidered in choosing the search method, The platform 
with the attached part, and the actuating mechanisms 
may have significant inertial moments which, in 
general, will not be identical for each of the n de- 
grees of freedom, 


Backlash and friction will be present in the 
Fig. 4. motion-transmitting mechanisms; inertial and resist- 
ance moments may depend greatly on the deflection 
angles of the actuating mechanisms, Random interferences may affect the search system; these are mainly caused 
by surface unevenness of the mounted blank. 


The known methods of automatic search in terms of many variables can be divided [2] into blind search 
methods (scanning, random search), and methods using an analysis of intermediate results (gradient and steepest 
slope, cyclic search of particular extremes (3]), The scanning and random search methods are associated with a 
long search time, and thus can be recommended only for the case when there are several extremes of the quantity 
being optimized, Q(x,,...» Xn), and organized search methods, which assume that Q has only one extreme, are 
inapplicable, The gradient and steepest slope methods require that not only the sign, but the magnitude of 
8Q/8 x; be determined beforehand, Because of the peculiarities of the systems mentioned, large error is incurred 
in determining the magnitudes of the partial derivatives, and thus these methods lose their basic advantages, 
Moreover, the latter require additional time for trial probes, and comparatively complex apparatus, The gradient 
method, in which each of n inputs is excited simultaneously by sinusoidal signals of different frequency [4], re~- 
quires that the components of the output signal Q, caused by each of the inputs, be determined, In general, a 
platform with six degrees of freedom means that the response Q to six input signals must be determined. Due to 
nonlinearities in the system, which cause higher-order harmonics, and because the inertial parameters of the ele- 


ments are not constant, analysis of the output signal becomes practically impossible, 
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Fig. 5. Fig. 6, 


The method of cyclic search of partial extremes in terms of each of n inputs is the most suitable method 
for determining optimum blank positioning (Gauss — Zaidel method), This method, in combination with the 
method given in [6], allows one to combine trial and working movements. For each of the variables in this case, 
only the first step can be made in an incorrect direction, and when the direction of movement is changed, back- 
lash must be taken up only once, Since the method requires that only the sign of 8Q/ 0x; be found, the peculiari- 
ties of the search system do not have as strong an effect as in the methods which have been examined, The given 
method permits one to accomplish search with very simple control apparatus, It will be shown below, however, 
that this method, without additional perfection, is not useful because of large error and long search time. 


Let us evaluate methods of searching for partial extremes in terms of each of n inputs from the viewpoint 
of the indicated criteria, We will examine the following methods of searching for the extreme of a function of 
one variable; 


1) self-oscillation method (5); 
2) with memory of the peak values of the extreme index [5]; 


3) with continuous change of the inputs, and measurement of the output increments at equal time inter- 
vals [6]; 


4) step method with constant time intervals between steps [5]; 
5) step method with variable time intervals between steps [6]; 
6) step method with constant step amplitude (5); 
1) step method with variable step amplitude [7]; 


8) extrapolation methods [8]. 


Methods 1, 7, and 8 require that the magnitude of 0Q/0 x; be determined, and, as a result, have a number 
of deficiencies which were mentioned above, Method 2 is not useful because of the effect of interference; indi- 
vidual peaks and valleys may be taken for peak values when the blank is moved, False change of the search di- 
rection results, and this leads to an impermissible length of search time. The advantages of method 5 cannot be 
used in the given system, since large time intervals between steps near the extreme lengthen the search process, 
Method 3 is the most suitable for the problem being solved, and also method 4 in combination with ~ethod 6, 
We note that the optimum action form is comparatively complex in methods 4 and 6, due to the presence of in- 
ertial elements and restrictions on the coordinates, The problem of determining it will be discussed in the next 
article, The methods given permit one to solve problems with comparatively simple control apparatus, 


SKE 


6, Effect of the Static Characteristic of the Search Object 


on Accuracy and Search Time 


The search method which involves finding partial extremes cyclically depends essentially on the distribution 
of the surface on which the extreme points lie relative to the coordinate axes. Thus, when the minimum of the 
function 


(12) 


is found by cyclic search along axes x, and x2, the system moves to the desired point O after only one partial ex- 
treme (O,) is found (Fig. 4a), If the symmetry axis of an elliptical paraboloid (12) is rotated with respect to the 
search axes, the search is substantially lengthened (Fig, 4b), The system proceeds consecutively from point A 
through the partial extreme points 1, 2,3,... 


Self-oscillations near a point lying to one side of the extreme are possible in such a system at the end of 
the search process, On the x,, xX, plane, let the Q = const lines have the shape of “elongated” ellipses, whose sym- 
metry axes do not coincide with the search axes (Fig. 5), Let us examine the fourth quadrant of the plane, The 
geometric locus of points at which the tangent to the ellipse Q = const is vertical is the straight line ON; the geo- 
metric locus of points at which the tangent to the ellipse Q = const is horizontal is the line OP, Therefore, in 
moving along a vertical line (for example M,, Ms), the partial minimum of Q is reached at the intersection of 
this vertical line with ON, and in moving along a horizontal line (for example Mg, My), the partial minimum of 
Q is reached at the intersection of this horizontal line with OP, 


Let us assume that the system moves in steps from M, as a starting point, At first, x, changes, and the sys- 
tem takes a step to point Mg, for example, Since the increment AQ > 0 in this case, a reversal occurs, and the 
next step is taken in the opposite direction; the system proceeds to point Ms, which coincides with M,. The sys- 
tem simultaneously remembers that AQ has been greater than zero once, If the same sign of AQ occurs in another 
step, the cycle of changing x, is completed, and a cycle of changing x, is begun, The system takes the next 
step from point Mg in the Ms, Mg, direction, Since AQ < 0 here, subsequent movement is in the same direction, 
In the next step, the system goes through point Ms. However, line ON is intersected in this case, and the system 
passes through a partial minimum, Let AQ > 0 in this case, Then searching is done along axis x,, moving at first 
along lines Mg, Mg, then reversing and moving along lines Mg, M7, and then M;, Mg. A partial minimum is passed 
through, since line OP is intersected, Let AQ > 0 in the next step, Then searching is again done along axis x9; 
movement is along lines Mg, Mg, then Mg, My, and My, Myy. The search is once more changed to the x, axis, 
moving along lines My, My, Mgg, Myg, etc. It can be seen from Fig, 5 that the mapping point of the system 
moves away from the true minimum point, and may move away a considerable amount, as will be shown below, 
One must use.a method which includes an alternation in the direction of the first step in the cycle so that the 
“staircase” My, Ms, Mg, My, « « - Will not lead the mapping point far from the minimum O, 


In Fig. 5, the first steps in the cycles for the variable x, are always taken in a positive direction — these are 
steps My, Mg; Mg, Mg, etc. The first steps in the cycles for the variable x, are always made in a negative direc- 
tion — these are steps Ms, Mg; Myy, Mg, etc. However, this is not obligatory, The first step in the first cycle of 
the x, variable may be made in a positive direction, the first step in the second cycle in a negative direction, 
first step in the third cycle positive again, etc, The alternation of signs in the first steps is analogous for x,. One 
can then see, by plotting a graph like that of Fig. 5, that the mapping point of the system will move in the vicin- 
ity of the minimum O, and will not get far from there, The generalization of this principle to any number of 
variables is evident, We will show that the deviation of the point around which self-oscillations occur from the 
extreme point may become significant if first-step sign alternation is not used, 

The lines Q = const of an elliptical paraboloid (12) are shown in Fig. 6, Let searching be started along the 
X, variable from point X- Search direction change, according to the principle mentioned above, does not occur 


at the partial extremes (denoted by circles in Fig, 6), but at the points Ag(Xs, Ag(Xg, Ra), 
Let pointsA,, Az, As, A, deviate from the partial extreme points by an identical amount equal toc, Then 


Is = +c, = 2 
1 as) 
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Fig. 7. 
In the region for which 


Ax, Ax, ={z,|—|z,| <0, 


(14) 


the system moves to the extreme point in the search process, When Ax, > 0 and Ax, > 0, the system moves away 


‘from the extreme. 


The equation of an elliptical paraboloid with its symmetry axes rotated by an angle is written in the form 


Q a> (x, cos a — 2 sin a)* + (x, sin a + 22 cos a)? : 


a b 
From the equation 
_ 
where x_ = X%_ = Xq. Xp = and the equation 
Fa, = % 


where x, = %, x, = % . we find the coordinates of points A,, Az. As, and Aq: 


= +0 = + 


where 
(b—a) tga 
A= b+atgta’ 


(b— a) tga 
p= 


(15) 


(16) 


(17) 


(18) 


(19) 


16) 


(17) 


18) 


(19) 


Using (13), (18), and (19), we find 


Z = ABz, + Ac +c, (20) 
= + AB% + ABc + Be+c. (21) 


We find the coordinates of the points around which self-oscillations occur (x4, Xg9) from the equations 


Az, = OandAz, = 0, (22) 
e(A+1) 
- B+1 


Leta = 2/4, b/a=10, Then = = 5.5 c, 
The given system will self-oscillate around the point with the coordinates which have been found, 


If c has a negative sign, the system, according to (23) and (24), will self-oscillate around the point with 
coordinates (Fig. 7) 


= = — 5,00. 


Thus, the point located in searching for the extreme of a two-variable function will be at a considerable 
distance from the true value, This distance can be greatly reduced by two methods: 


1) by changing the optimizer structure so that the sign of the deviation of the system from partial extreme 
points is changed in a repetitive cycle (by alternating the direction of the first step along each of the variables); 


2) by putting a coordinate transformer in the system, whichchanges the direction of the search axes, 


Let us examine the first of these methods, In the method of alternating the first-step direction, equations 
(13) will take the form 


= 


om 
Self-oscillations will occur around the point with coordinates 

— (26) 

— (27) 


As in the case examined earlier, let a = 1/4,b/a=10, Then xj =1.1c, xf =1.1c. 
In this case, the coordinates of the point located will be near the extreme point. 


It follows from (26) and (27) that the search error can be greatly reduced by bringing & to zero or, in prac- 
tice, by diminishing it. This problem will be treated in detail in the next article, 
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Original article submitted June 2, 1960 


This article presents a method for solving sets of finite equations by means of analog computers, 
An auxiliary set of differential equations is constructed from the initial set of finite equations, 
and differs from them in that the stable rest points correspond to the roots of the initial set. In 
this paper, the particulars of this method are studied, and the conditions under which it is most 
effective are found, 


1. Introduction 


Modern mathematical indirect analog computers are basically intended for solving various dynamic pro- 
blems from the class of ordinary differential equations, both by the nature of their general structure and by 
their individual elements, 


This, however, is not the only type of problem which can be successfully solved by these devices, It is 
known that it is possible, in principle, to solve algebraic and transcendental equations, arising from the solution 
of certain differential equations under special conditions, by using such computers. 


The method of solving algebraic and transcendental equations by such means will henceforth be called 
the “differential equation” method, Historically, this method arose as a result of work on the elimination of in- 


stability which occurred in solving algebraic equations by the implicit function method on continuous action 
machines, 


According to this method, instead of using the initial equation 


f(x) =0 a) 
one uses either its equivalent form x = x + f(x), or an equation which is nearly like Eq. (1), 


k>1. (2) 


A block diagram of the problem in the form of Eq. (2) is shown in Fig. 1. It follows from Fig. 1 and Eq. (2) 
that the structural diagram is a closed-loop network, and that the dynamic properties of the elements are not con- 
sidered, In view of this, the solution x* = (a is a number), which tums Eq. (2) into an identity, may prove to 
be unstable, 


Fig. 1 Fig. 2 

Thus, the problem arose of how to ensure network stability for sufficiently large k. This problem was com- 
plicated by the fact that it is difficult to account for small inertias which reflect the dynamic properties of the 
network, and in practice, a differential equation was solved which did not correspond to the block diagram of 
Fig. 1. 


However, one remarkable property of dynamic systems has been known for a long time. If one large con- 
centrated inertia is artificially put into such a system with small inertias, the system will be stable in the major- 
ity of cases, if the differential equation formed without accounting for the small parameters has a stable solution, 
The latter statement evidently requires some explanation, In the case being considered, the artificial inertia is 
created by replacing the high-gain amplifier with an integrator (Fig. 2). In this case, the differential equation, 
neglecting small parameters, is 


=—f (2), @) 


where T > 0 is the integrator constant, Substituting T = t/T, we get 


F =10). 


The hypothesis of the suppression of small inertias can be formulated as follows: If Eq. (3") has a stable 
solution x* = ©, which is the root of Eq. (1), then a number T # 0 can be found for which the actual system shown 
in Fig. 2 will be stable at the point x = x*. This definition gives rise to the requirement that the rest point x* of 
Eq. (3") be stable, This equation is generally formed without considering this requirement, since its right half is 
arbitrary. 


The hypothesis may be untrue even if the latter condition is met. However, practice shows that if this con- 
dition is met, a number T for which the system is stable will exist for a broad class of problems, The method pre- 
sented for solving algebraic equations is the differential equation method, At the equilibrium point dx/dr = 0, 
Eq. (3*) turns into Eq. (1), and we get the solution x = x* at point a (Fig, 2), The derivation presented gives rise 
also to an inverse relation: the differential equation method is inapplicable in the form described here if the 
rest point x* of Eq. (3") is unstable, Practice shows that the majority of equations fall into this class, 


A new problem arises: How does one construct a differential equation which will have the same equilibrium 
points as Eq. (3), and have all these points asymptotically stable? 


Methods for constructing such equations are known in the literature [1], The following equation: 


dz af 


is written instead of Eq. (3). 


The right half of this equation is a derivative of x having the form A f*, and the equality itself shows that 
the form is minimized in terms of the gradient, No detailed study of this is given in (1]. 


This article is devoted to a more detailed study of this method, The peculiarities of the method are ex- 
plained, the region of its most suitable application is found, and stability in the large is examined, It will be 
shown that the gradient method, as it applies to the stated problem, is derived from Lyapunov's second method, 
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2. Statement and Solution of the Problem 


Let a set of algebraic or transcendental equations exist: 


=O (i a), (5) 


where Xs, Xg, + « « » Xp are the variables for which the given set must be solved, The functions fj are continuously 
differentiable over their entire arguments, at least in the bounded region D. In the vicinity of the solution 
(3, 


det A=det} -----.-- +0. 


The latter condition ensures that the rest point will be isolated. The determinant (6) may equal zero on 
individual manifolds, These cases, however, will be dealt with specially, 


_ We will use Lyapunov’s theorem about the asymptotic stability of the equilibrium point [2] to construct a 
set of differential equations having the properties we need, According to this theorem, the equilibrium point 
+ » %p) of the set of differential equations 


d: 


is asymptotically stable for sufficiently small values of y; (i = 1,..., 1) if a positive-defined function* 
V = V(yq, «++,» Yn) exists whose derivative with respect to t is negative-defined. 
The content of the theorem can be written as the inequality 
(8) 
i=1 


Let us form, by some means, a positive-defined function V of f;,...,. fn: 


V=V(fi,... 


’ Tn); 


t+ Sn + Yn) (i=1,..., a). (9) 
We note that it will also be positive-defined in terms of the arguments y;,.... Yn« 


The function goes to zero if f; = fz =... = fn, that is, at the point (x1,...,Xn) = y,..., %n), where 
xj = Oy + yg (i =1,..., 2), which is a root of the initial set of finite equations (5), 


It follows from inequality (8) that, if one forms a set of differential equations of the form 


dy ov 


‘on (10*) 


the conditions of the theorem will almost be met,**that is, a positive-defined function V exists whose derivative 
equals 


*A positive-defined (negative-defined) function V is a function which is positive (negative) everywhere except 
at the points y; = yz = . = Yn = 0, where V = 0, 
** Strictly speaking, the conditions of Lyapunov's other theorem — about nonasymptotic stability — will be met. 


| where 
or its equivalent 
dz; av 
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or at least is nonpositive: 


V<0. 


It follows from this that one cannot always guarantee asymptotic stability even in the small, However, if 
we use the fact that the coordinates of any physical system are always subject to the effect of small random forces, 
it can be shown that asymptotic stability will exist in the majority of cases. 


Indeed, let dV/ dt = 0 on the manifold M of dimensionality m < n, i.e., a dimensionality less than that of 


the space, Then a positive probability exists that the mapping point will be in space L,, formed from L by being 
taken from the manifold M, The inequality 


dv 
ar < 9 (11) 


holds strictly in the latter space, 


If the mapping point, affected by small perturbations, is thrown onto a trajectory moving away from mani- 
fold M, Lyapunov's theorem about asymptotic stability will be justified from there on, Let us clarify this fact by 


an example of a space where n = 2 (Fig. 3), In terms of Eq. (7), let dV/dt = 0 on the manifold V =, and, con- 
sequently, the set of equations 


reduce to the identity OV /0x; = 0 everywhere. on this manifold, 


Let us assume that the trajectory intersects V = C at the point (xf, xf). Because of fluctuations, the map- 
ping point can be in the €-vicinity of this point. But. inside V < C, V< 0, and the point, once thrown there, will 
move away to the equilibrium position, 


It is evident that none of these statements is justified if dV/dt = 0 on a manifold of dimensionality m = n, 
or over the whole of space L. 


Thus, the use of the hypothesis of small perturbations greatly weakens the conditions of the cited Lyapunov 
theorem, 


Before considering a concrete choice of the function V, it is necessary to examine some possible cases of 
system behavior for xj €D, where D, as mentioned earlier, is a bounded region, This examination is obligatory, 
since it is known [3] that the system behavior near the isolated special point does not fully determine the be- 
havior in a sufficiently large, though bounded, region. 


It is evident that the system will be asymptotically stable in the presence of any perturbations xjp €D, if 
the structure of V is the same everywhere as in the vicinity of the rest point, Near the rest point, the function V 
in i-dimensional space is a one-parameter family of closed surfaces enclosed in one another, 


¢ 
é A | 
J q 
y 
v=C 
Fig. 3 Fig. 4 : 
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It is easy to give an example where such a function is destroyed if r = })| xio| is sufficiently large. 
i 


Figure 4 shows a function V = V(x) > 0 for x # x*, and V(x*) = 0. Let us examine the stability of point A, We 
will locate the origin of the coordinates at this point. 
We will form Vy = V(x, + x)— V(x,), which is greater than zero, and V(x, + 0) = 0, Let us examine its de- 
rivative Vj, in terms of the motion equations x = — 94/2 x; 
‘av 


where x = |€] and Vj, (x4) = 0. Asymptotic stability of point x, will consequently be observed in the presence of 
perturbations not exceeding those mentioned, Points B and C are evidently unstable in the presence of fluctua- 
tions, Similar arguments can be presented for any finite n, As we see, points A, B, and C are stationary points 
which are defined by the set of equations 6V/ Ox; = 0 (i=1,...,1). Consequently, one must attempt to con- 
struct the function so that it is geometrically a “smooth cup" in (n + 1)-dimensional space, with a minimum at 

the point under investigation for stability. It is true that this function is ideal, A function can be constructed 
which has “safe* extremal points like point C (inflection points, saddles, etc,), For this is needed partial derivative 
8v/ 8x; #0 for any xj, except at points corresponding to the solution of Eq. (5), 


The discussions on the construction of the function V presented above can be used in actual cases of equa- 
tions (5), and are also needed for further derivations, 


We note that one can choose a positive-defined* quadratic form 
n,n 
V= > jf if aij = const 
i,j 


for the function V., The canonical form, in particular, may be used: 


(12) 
i 


It is easy to prove that if Aj = 1(i=1,..., 1), the analog set of equations 


(14) 


where dx/dt and f are column matrices, A‘ is a matrix which is a transponent of the Jacobian matrix A. It fol- 
lows from Eq, (14) that condition (6) is met at the solution, and all false singularities are on the manifold 


det A=0, 
False solutions need not be feared if the analog is universal, since they can be isolated, If the solution is 
n 
false, >} | /4| £0. Im the case of a specialized analog, when the path of the root , as the left half of (5) changes 


with time,is important, an “ideal* function V must be constructed by some means or other, or else the root may 
go off on a false path, 


Let us examine some specific properties of Eq. (13), Equating the right halves of (13) and (14), we get 


*A positive-defined quadratic form is a quadratic form whose matrix is positive, that is, all the main minors of 
the matrix are greater than zero, . 


— (i=1,..., 2) (13) 
t 
can be written in the form 
dx 
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ve 


where 0V/0x is a column matrix of partial derivatives, Let us assume that 


det A #0 


in the region bounded by V =C (except at the rest point, where the determinant can go to zero), Then it can be 
demonstrated that the column matrix 


A'f +0. 


Indeed, the matrix A‘ is nondegenerate if only one of the fj is not zero in the vicinity of the rest point. 
Consequently, in view of identity (15), only one function 8V /0x; is nonzero everywhere in this region, 


Remembering that 


we get the final result v< 0, and this means that the conditions of the asymptotic stability theorem are strictly 


met, The point is asymptotically stable in the presence of perturbations not exceeding the limits of the region 
being examined, 


Let us take an example, Let there be a set of transcendental equations of the form 


fy a, = const, 
fo = e+ a,=0, == const. 


It can be shown by graphical construction that this set of equations has a real solution for certain values of 
a, and a, (for example, a, = a, = 1), We will explain stability in the small, To do this, we find the roots of 


Of, 


af afs 
=iti 


The roots have positive real parts, so a set of differential equations of the form 


has an unstable equilibrium point; We will construct a new set of differential equations by using equality (13), 
We will use V = $(f% + #3) as the function V, Then 


. a a 
Substituting the expressions for f;, we finally get 


= — | (ty — + 1) + + ag) 
= — [— — + + (€* + + 


In the example being considered, det A # 0 for any values of (x;, x,). Indeed, 


dot A = 1 + ett 0 


av 
=4), (15) 
| 
. 
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and, consequently, the rest point, being a solution of the transcendental set, is asymptotically stable in the large ,* 
that is, for any initial conditions, 


Equation (14) comes as a result of Lyapunov's theorem, This theorem gives the sufficient conditions for the 
existence of stability, and thus it is desirable to try to simplify this set of equations, Practice shows that one of 
the simplifications is to make all the elements of matrix A zero except along diagonals; this often gives positive 
results, The analog system in this case assumes the form xj = —0 fj/ Ox; (i=1,...,m). This is significantly 
simpler than (13), 


3. Finding the Roots of Analytic Functions by the Gradient Method 


The statements presented in the preceding section on the properties of an “ideal” function V permit one to 
state the following problem, Let a function V of f; (i= 1,..., 1) be given by formula (12) (Aj =1;i=1,... 
«++,M), It is necessary to find a class of function which, when substituted in this formula, will give an “ideal” 
function V. It happens that one of these classes is that consisting of the real and imaginary parts of any analytic 
function of one complex variable. Let a real analytic function f(p) of complex argument p=&+iB be given. 
The problem is to find its roots % +iBy(k=1,..., q) by the gradient method. To do this, we form the function 
+ iBy, 


V = + (a, 8) + N*(a, BY], (as) 
where M(a 8 ) and N(a ,6 ) are functions corresponding to the representation of the initial function f(p) in the form 
{(% + iB) = M (a, B) + iN (a, 8B). 
The function f goes to zero, corresponding to the root &; +i8;, when, and only when 


M (a, 8) 0, 


N (a, B) = 0, (17) 


Let us solve (17) (which is generally nonlinear) by the gradient method, We will first minimize (16), It is 
evident that it will be zero when, and only when, equality (17) is satisfied. Let us examine the set of differential 
equations 

da aM an dB aM an 


dt da a ap (18) 


which are formed from (10). 


It follows from Liouville*s theorem [4] that V is an infinitely large function (V + @ for |a| + |8|— o), 
since every analytic function possesses this property, providing it is not a constant and is analytic everywhere. 


This fact permits one to state that, for any initial conditions, the mapping point goes to one of the mini- 
mums of V during a course of time, and goes to one of the roots of this function if there are no false minimums; 
that is, minimums which do not correspond to roots of f(p). 


If there are no false roots, V is shaped like a series of one-sided cups with minimums at roots, intercon- 
nected by saddles, The entire analytic region in the (a ,B ) plane is divided into a series of contiguous subspaces 
connected at their roots, Let us proceed with the proof itself, Suppose that the function V has no minimums or 
maximums at points defined by the inequality f(p) # 0, and that a two-d#mensional continuum of stationary 
points of the other type does not exist, We note that the function V is one-half the square of the modulus of f(p). 
It follows from the principle of the maximum modulus of a function of a complex variable that V has no maxi- 
mums or minimums at any point in the analytic region except at roots, where | f(p) | = 0. We will now show 
that the analytic function does not have a continuum of points at which 


*In judging the stability in the large , one can see that the function V is “infinitely large" [3]. A sufficient con- 
dition of an infinitely large function is that V - o for E |e |— =. This condition is met in the example pre- 
sented, 
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Let condition (19) be met in some region D. It can then be proved that V is constant on this continuum 
(see Appendix), It is known [4] that, if the modulus of an analytic function is constant somewhere in a region G 
belonging to the analytic region, then V = const over the entire analytic region, We have a function V # const, 


In this case, the derivative ie --|(%)+(%)| goes to zero on a manifold of dimensionality less than 


that of the space, which is what we wanted to prove, Thus, it has been established that the solution will con- 
verge to one of the roots of the function f ( + i8) if there are no perturbations in © and 6 exceeding the limits 
of the analytic region, It follows from this that, in using the gradient method, the problem of finding the roots 
of analytic functions reduces to one of finding all contiguous regions; these completely fill the (a ,6) plane. 


As is known, integral functions, polynomials in particular, are related to analytic functions, Consequently, 
the well-known problem of finding the roots of a polynomial can be solved by the gradient method, 


The case of finding the complex roots of one equation has been examined above, It would appear that a 
set of differential equations will possess the property of not having false rest points even in the multidimensional 


case, if we take V= Bn) On, Br, Bn)] a8 the Lyapunov 
i 


function, Here, Mj and Nj are the real and imaginary parts of the ith equation of the set fj(py,.... Pp) = 9, 
respectively, It is easy to prove, however, that such a function is not the modulus of an analytic function of many 
complex variables, and thus one cannot guarantee that the set 


at da,’ R, =1,..., a) 


will not have false rest points, We will mention here one method which may give positive results, Let us con- 
sider the set of differential equations 


For any ith subsystem of equations, the variables «,.,8, (k # i) are independent parameters, which deform 
the function V in time, The movement of the mapping point in the plane of subspace (4,64) is a function of 
the gradient. However, if the rate of deformation is significant, the movement may cause the rank of V = Cj to 
increase, i.e., instability will be observed, 


4. The Effect of Bounding the Linearity Limits of the Analog Elements 


on the Stability of the Whole System 
Let us examine the block diagram of a network for solving a set of two equations® (Fig, 5), It consists of 
two functional converters @, = — 5” and@, = — <— and two integrators. Two loops I and II are clearly 
2 


shown in the network, corresponding to the differential equations X, = @, and x, = @,. Initial conditions are in- 
troduced in the usual way — capacitances C, and Cz are charged to Xy and X99 (with switches K, and K, open), 
After closing K, and K,, two modes must be distinguished: 


1) integrators I, and I, operate in the linear region when t > 0; 
2) the integrators operate outside the linear region part of the time, 


The second mode is related to cases when the absolute magnitudes of the coordinates | x |, for t > 0, at 
certain times assume values greater than the initial ones, The functional converters may either operate linearly, 
or may depart from linearity, in both cases, 


*If n > 2, the network remains similar, except that the number of closed loops increases in proportion to n, 


eV. 
Vi= Mi4+M,, 
av 
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ae Fig. 6. 
Let us investigate both cases from the standpoint of the stability of the rest point defined by the equations 
,=0,,=0, 


We will examine the first mode, Let the linear zones of the functional converters be max | ®;|=A;, Aj > 0 
(i=1,2), It is also known that if |®,| <A, and |,| <A2, the process converges to the roots (x,*, xf ) of the equa- 
tions 6, = 0, 6, = 0. We are only interested in the behavior of the system in the nonlinear region, Let us examine 
the function* 


where f, and f, are the initial functions of the unconverted system f,; = 0, f, = 0, and its time derivative is 


= 95,71 +5 % 


But Acsign oY — 0 and sign = —1 when 
Ox; 


V = Ay sign Ay sign 2 me) <0. 


Thus, the first mode is safe, and convergence of the process to the roots of the initial system is ensured. 
We will examine the second mode, assuming that the nonlinear functions £, (€ ;) are smooth and approach the 
saturation boundary |C| asymptotically (Fig. 6). We will write €, in the form 


The function M(E ,) is positive everywhere, and approaches the axis of asymptotically as We 
will assume that a block with a transfer function M(#;) is connected in series with the functional converter 0; . 
Let us write the differential equation of the integrator, taking into account the nonlinearity of the differential 
coefficient k, As is known, the transfer function of the integrator has the form 


k d 
k>0, 


We will assume that the amplification factor k is a function of the output coordinate x, in the same man- 
ner as M(€ 4). Consequently, the integrator equation is 


dz x 
*It is assumed that V has the "ideal" form. 
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where y is the input and T is the integrator constant, The analog set of equations is then written in the form 
M (i = 1, 2), T =14., (20) 


Considering (20), the derivative of V with respect to t can be variable in its sign: 


2 
wy 

For large xj, i.e., for significant departure from the linear zone, the second term may be large, and may 
determine the sign of the derivative V. This kind of mode is dangerous — false rest points may appear (for ex- 
ample in the saturated state), or self-oscillations may arise, When this mode appears, the scale of x; must be 
altered in such a way as to decrease the range of change, 


One can draw the following conclusions on the basis of the material presented; 


1. In solving a set of finite equations with an analog computer, it is desirable to construct the initial dif- 
ferential equations by means of the gradient method, This gives a logical approach to the composition of the 
differential equations, and ensures convergence to the rest point, which is the solution of the set of finite equa- 
tions, in the majority of cases, 


2. The rest points of the set of differential equations must be asymptotically stable, It has been shown in 
this paper that, due to the positive effect of small perturbations, the gradient method ensures asymptotically 
stable rest points in the majority of cases, 


3. A sufficient condition of asymptotic stability of the equilibrium point when V is expressed in the form 


n 
V=+ > , and when the perturbations do not exceed the limits of a contiguous region, is that the deter- 
a 


minant of the Jacobian matrix be nonzero throughout this region. 


4, Asymptotically stable rest points of the set of differential equations may include some which are not 
roots of the initial set of finite equations, 


5. It has been proven in this paper that, if the left half of a given finite equation is an analytic function, 
a set of equations whose independent variable may be complex will not have asymptotically stable false rest points, 


6. If integrators and functional converters operate nonlinearly, two operating modes must be distinguished: 
a safe mode, in which only the functional converters are permitted to operate nonlinearly, and a dangerous mode 
in which the integrators operate nonlinearly, Stability may be lost in the latter case, 
APPENDIX 


If a differentiable function V(x,y) meets the conditions 8V/ x = 0, PV /dy = 0 at all points in a simply 
connected region D, then V = const in this region, Let us examine two different points (x,y), (x",y"), and AV = 
= V(x",y")— V(x,y). The points (x,y), (x*,y*) are so chosen that the triangle (x,y), (x,y"), (x",y") lies entirely in D. 
Let us assume the contrary, that is, that V + const, Then AV #0, We will form the incrementAV in the fol- 
lowing way: 


Using Lagrange's theorem, we get 
where x < < x" and y< < y*, The points (€ ,y"), x.) €D and, consequently, 


az = oy 


| 
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It follows from (21) that (x*- x) = (y"— y) = 0, which contradicts the initial assumption that two different 
points were chosen, Thus, the theorem is proved by reductio ad absurdum, 
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This article presents a method for choosing the optimum structure of digital analog computers 
on the basis of a special factor — the figure of merit. Comparable structures are described by a 
series of numbers called the structure parameters, The problem is to find some group of param- 
eters which give the extreme value of the figure of merit, A concrete example of a design of a 
digital analog computer is given. 


Specialized computers intended for operation in automatic control systems must satisfy a number of spe- 
cific requirements on accuracy, speed of response, reliability, size, operating convenience, etc, To a great de- 
gree, satisfying these requirements depends on the rational choice of a structure for the computer. Therefore, 

a comparative evaluation of computer structures and the problem of an optimum structure in one or another man- 
ner are of direct interest in rationally constructing control apparatus, 


A special class of computers, called “digital analog" computers, deserves attention in connection with the 
solution of control problems [1,2]. Digital analog computers are distinguished by a number of advantages in ap- 
paratus and use, On the other hand, the majority of computers of this type described in the literature, which are 
the series type, have a slow response, Attempts to raise the response speed by using parallel and series- parallel 
structures lead to a rapid growth in the number of elements which, in turn, makes the apparatus more expensive 
and requires additional measures to ensure reliability, These contradictions make the problem of choosing the 
optimum structure of digital analog computers especially vexing. 


This article is devoted to an examination of this problem. However, it does not pretend to be a complete 
investigation of the subject, and is only one of the first steps in this direction, 


Statement of the Problem 


Speed of response, accuracy, and reliability of computers depend primarily on the characteristics of the 
elements used in them, However, in the majority of cases, the computer designer is limited by the elements 
readily available, chosen on an economic basis, or representing the best made available by the state of the art, 
Under such conditions, the required characteristics of the computer may be obtained only by using a number of 
elements, and by suitably interconnecting them, that is, by changing the computer structure, 


Henceforth, we will mean a system composed of standard elements working as a single computer when we 
speak of the structure of a computer, Standard elements, in turn, mean those which are comparable in complexity, 


construction, and function; these are monostable multivibrator summators, memory cells with aone-bit capacity, 
logic control elements, etc, 


By varying the number of elements in the structure, it is possible to get a relative enhancement of one 
computer characteristic at the cost of worsening others. Thus, for example, increasing the number of parallel- 
acting computer channels*® ordinarily increases the response speed of accuracy, but simultaneously lowers the 
reliability and makes the computer more difficult to use, Therefore, it is natural to compare different structures 
on the basis of an over-all examination of the characteristics which are basic to the given type of computer, 
thereby determining the relative efficiency of the structures, If the structure efficiency is evaluated in some 
special way by a quantity — the figure of merit — equivalent structures will have the same figure of merit, 


On the basis of these considerations, the problem of choosing the optimum structure of a digital analog 
computer is examined below, By optimum structure, we mean one which satisfies the given technical require- 
ments, and which is characterized by an extreme value of the figure of merit, 


In the first part of this paper we present a choice, and the basis for it, of the figure of merit for a structure, 
In the second part, a method for choosing the optimum structure, using the figure of merit which will have been 
derived, is given, The third part deals with an examination of a concrete example of choosing the optimum 
structure, Henceforth, it will be assumed that an algorithm for the structure is given, 


Figure of Merit 


The basic task of a computer in an automatic control system is to process information fed to it from cer- 
tain parts of the system, The results are transmitted to other parts of the system and are used in the over-all 
control process, From this standpoint, the computer can be considered as a communication channel, and its 
quality or efficiency of structure can be evaluated as the traffic capacity of such a channel, Of course, we con- 
sider that the traffic capacity is, to some degree, an idealized characteristic, and that the real amount of in- 
formation transmitted by the channel per unit time depends on a number of factors, primarily the type of pro- 
blem solved by the machine, the computer program chosen, etc, These factors, in turn, depend on the way the 


computer is used, as well as how rigidly the computer links the input and output signals, independently of their 
information content, 


The traffic capacity increases as the number of parallel computer channels is increased and, consequently, 
cannot fully characterize the quality of the computer. Considering this, computer quality or efficiency should 
be evaluated by means of its traffic capacity referred to some quantity which characterizes the amount of ap- 
paratus used to get such a capacity or, in other words, by the apparatus complexity of the computer. 


The betrer of two structures having the same traffic capacity will be the one which has the fewer number 
of elements, other conditions being equal, However, different types of standard elements are not completely 
equivalent from the standpoint of reliability, cost, convenience of use, etc, Hence, the apparatus compiexity 
of a computer is more correctly evaluated by an expression of the form 


= + +.. AyMy, (1) 


where nj is the number of standard elements of the ith type used in the computer, and a; is a coefficient selected 
on the basis of experience with analogous computers, or other available data, 


The a; coefficients can be made proportional to the probabilities of failure of each of the standard element 
types. If one then assumes that redundancy is not introduced into the structure to increase reliability, the quan- 
tity A will be proportional to computer failure and, consequently, will be inversely proportional to its reliability. 
In the literature on reliability [3,4], the X computed under similar conditions has been called the complexity co- 
efficient, and relation (1) has been called the complexity relation, 


The material considered permits one to express the figure of merit of a computer in the form 
where c is the traffic capacity of a computer which has a given structure, 


*We mean here an increase in the number of parallel channels not intended to introduce redundance, 


| 
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If, in choosing the figure of merit, not only the traffic capacity of a structure, but its reliability is con- 
sidered as well, and the possibility of raising reliability by introducing redundancy is taken into account, the 
figure of merit will certainly be different. The situation becomes more complicated if the cost of building and 
operating the machine is considered along with reliability. However, it can be assumed that both reliability and 
cost, if given types of standard elements are used, will be determined basically by the number of these elements, 
and expression (2) can be taken as the simplest of a number of analogous evaluations. 


Method for Choosing the Optimum Structure 


The method given in this article for choosing the optimum structure is as follows, After the figure of 
merit is gotten, some initial structure is composed which satisfies the given algorithm and is used as a first ap- 
proximation, A number of quantities called structure parameters are introduced on the basis of an analysis of 
this structure, One must keep in mind that the parameters chosen must describe the structure accurately enough 


so that it may be distinguished from other possible structures which satisfy the given algorithm, That this may 
be done will be shown below in a concrete example. 


After composing the initial structure and choosing its parameters, the basic functional dependencies relat- 
ing the figure of merit to the structure parameters are derived, The possible changes which can be made in each 
of these parameters are determined by analyzing the technical requirements on the computer, Then the problem 


of choosing the optimum structure can be formulated in the following manner: find the maximum of some func- 
tion 


= G(F) (3) 
81 (F) > m, 


(4) 


where 7 is the figure of merit, F = ||F,,...,Fj{,...|l are the structure parameters, and mj are the boundary 
conditions, 


Such a problem may be solved by one of the existing methods, such as the gradient method [5]. In connec- 
tion with this, it should be mentioned that conditions (4) determine the region of change permitted in the struc- 
ture parameters, atid functions (3) can reach a maximum both inside this region and on its boundary, Under cer- 
tain conditions, the occurrence of a maximum on the boundary is sufficient cause for reexamining the initial as- 
sumptions, or for altering the technical requirements. 


The group of parameters obtained is adjusted on the basis of a number of additional considerations which 
either could not be expressed in (4), or were neglected in the basic solution in order to simplify the problem. 


An Example of Optimum Structure Choice 


Sor planing digital analog anslyser consisting of M independent blocks ~ 
integrators ~ which performed the trapezoidal integration approximation 


An) An, (5) 
i=1 


where z is the approximate value of (u(aae, yi is the value of y(x) at the end of the ith division of the x 


Xe 
axis, Ay; is the increment of y at the ith division of the x axis, and Ax; is the size of this division, After com- 
puter numbers [1] have been introduced, expression (5) is converted to a set of recurrent functions: 


= Yi-4 + + 
Ri+ Aa; (ti) = Via, (ti), 


g3(F) = ms, = 
| 78 


1 for Yi ax (4) > A, 
a, (ti) = 0 for 0< Rint (ti) <A, 
—1 for Y; ax (ti) << 0, 


5 ay (ti) (6) 
(in this case, A > R > 0), where Y is a number stored in the Y registers of the nopieeatts Y® is a number trans- 
mitted to the R register of an integrator; R is a number stored in the R register; a,(t;) is an input signal represent- 


ing increment of x; ay(t;) is an input signal representing a y increment; a7(t;) is the integrator output, represent- 
ing a z increment, and A is the maximum of Y and R (see-[2]). 


The outputs of each integrator occur as delta- modulation pulses with three possible values (—1, 0,1), The 
expression for the traffic capacity of such a channel was derived by V. M. Baikovskii. It is: 


C= (1 2cos bits/sec, (7) 


where A is the number of quantum steps in the signal, which equals, in this case, 


= 


(8) 


where N is the total number of stores in the y and R registers (see below) and 5, is the time between two 
successive outputs or the duration of one iteration, 


If the number A is sufficiently large, one can write, approximately, 


c= log, 3 bits/sec (9) 
t 


and the final form of the expression for the figure of merit of a digital analog analyzer consisting of M independ- 
ent integrators is 
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A block diagram of an initial structure which satisfies algorithm (6) is shown in Fig. 1. The numbers Y 
and R, increments a,(t;), and the program which determines the order of interconnection 
MB AD of the individual blocks are stored in memory block MB, The operations required by 
formulas (6) are performed in adder AD, namely, 3M successive additions at each step, 
The control apparatus C connects the blocks and sets the scale of the respective variables, 


Cc The operations of computing formulas (6) are always done in the same sequence, 
Therefore, address to the memory will always occur in a definite sequence, This elimi- 
Fig. 1 nates the need for an address system, and makes a delay line or shift register the best 


memory to use, The structure of the memory, under these conditions, can be described 
by the following set of parameters: the capacity F, of the part of the memory which 


stores Y and R, and the capacity F, of the part of the memory which stores the increments az(tj), the program, 
and the scale factors, 


In the simplest case, adder AD may be a single monostable multivibrator binary summator and a delay 
line which stores the transfer signal. In the general case, the adder will consist of a number of summators con- 
nected in parallel. It is also possible to have the three additions of (6) done by three independent summators 
connected in series. Starting from this basis, the adder structure may be described by two parameters: the num- 
ber of parallel channels Fs, and the number of summators connected in series in each channel, Fy. The para- 


meter Fs is also related to the memory, since the number of parallel channels in the memory must correspond to 
the number of parallel channels in adder AD? 


*The case may certainly arise where AD and MB have a different number of parallel channels, and appropriate 
converting elements are connected between them, However, this variation is clearly not rational for such a 
simple machine as a digital analog computer, and will not be considered here, 


= 


The control apparatus of the digital analog computer 
y eee: Aa is simple, and can be physically realized in a number of ways, 
=: In terms of the given analysis, it can be characterized by a 
single structure parameter, namely, the number of elements 
A used, Fg. We also need to know the values of the a; coeffi- 
cients in the expression for the complexity coefficient (1), 
SQ and also data on the speed of response, so that the time for 
x one iteration can be found; i.e., in essence, the character- 
TINE, <1 istics of the elements must be introduced, At present, such 
NA=3 complete data on all the different types of elements do not 
exist, Therefore, we will only consider certain basic types of 
elements here, and rough estimates of their characteristics 
will be used, 


+4 


i 


| 
{ 


AY 


£88 


Magnetic drums, delay lines (magnetostrictive or elec- 
tromagnetic), diode-ferrite shift registers, ferrite core matrix 
Fig. 2. memories, and transistor-capacitor memory cells will be con- 
sidered as possible memory elements, Only transistor-diode 
elements will be considered for logic circuit elements, in particular the control apparatus and adder, The unit 
for the number of elements is assumed to be an active element consisting of a logic cell which converts or re- 
members one bit of information for one cycle, Appropriate factors are introduced when other elements are 


examined. Thus, for example, 10 active elements per single information storage channel are assumed for a mag- 
netic drum, 


~ 
A 
a 


66 100 200 1200 


The response speed of the elements is evaluated in terms of the maximum permissible cyclic frequency 
for the memory elements, and the delay introduced for the logic elements, Thus, the maximum cyclic frequency 
is set at 100 kc for magnetic drums, magnetostrictive delay lines, diode-ferrite shift registers, and matrix memo- 
ries, In the case of electromagnetic delay lines, the maximum cyclic frequency is determined by the auxiliary 
electronic apparatus, and is set at 1-3 Mc, The maximum cyclic frequency is 1 Mc for transistor-diode memo- 
ries, The delay introduced by a transistor-diode logic element is estimated to be 0,1 [lsec. 


In finding the dependence of the figure of merit on the structure parameters, we also take into account 
the fact that when F; = M and F; = NM, where N = 2 + log, A — the full number of binary steps in X, Y, and R, 
including sign and additional minor rounding-off operations — the structure of the digital analog computer can 
be completely represented as an assemblage of individual functional blocks which are autonomous from the stand- 
point of physical realization. In this case, the individual blocks can be interconnected by busses, as in ordinary 
computer practice, and scales can be set manually by means of switches or similar devices. We can therefore 
assume that when Fs; has the values mentioned, F, and Fs go to zero, Finally, in evaluating the response speed, 
we take into account the fact that, if the number of parallel additions does not exceed N, the summators operating in 
parallel must be in a cyclic transfer loop, 


Figure 2 shows the figure of merit as a function of F, and F, for the cases when a magnetic drum or transist~ 
or-capacitor shift register is used as the memory. In calculating, the most frequently encountered values, M = 
= 60 and N = 20, were used, It can be seen from the graph that a pure parallel structure is better when shift 
registers are used, while a series-parallel structure with the number of channels equalling the number of inte- 
grators is preferable if a magnetic drum is used, The series-parallel structure is especially convenient if the re- 
sponse speed is not critical, 


A block diagram of a digital integrator, designed by means of the method described, is shown in Fig. 3. 
The integrator consists of N total binary summators £4-2y, in a cyclic transfer loop, and N shift registers P,-P}y 
with a capacity of two bits, One step of the computations in (6) is performed in eight cycles corresponding to 
the times when the shift pulses are fed to registers P. In the first cycle, the operation 


1 
is performed, the number Yj - ; is fed from the registers to the summator inputs A, and the number ay(t;) is 


shifted one step to the right, and is fed into the B input, Information is transmitted inside the registers in the 
second cycle, In the third cycle, the operation 
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R, + Aa, (ti) = Rin + Yias (4) 


is performed. 


The number Rj .; is fed from the registers to the summator inputs A. The number Yj enters the sum- 
mator inputs B through gates Ky~Ky if ay(tj) = 1, and through inverters Iy-Ijy and gates Kj-Kjy if a,(u) = ~1. 
If ax(tj) = 0, gates K and K* are closed, and R,-; is added to zero, The resulting signal az(tj) goes to the logic 
block, In the fourth cycle, information is transmitted inside the registers. In the fifth cycle, the operation 


Yi=Vi + 
is performed, 


— 
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Pp, 
CBA 4 
Input | 4@r(ti) 
Input block 
Fig. 3. 


The network operates in the same manner as in the first cycle. In the sixth, seventh, and eighth cycles, in- 
formation is transmitted inside the registers. Moreover, az(tj) = +1 signals are emitted from the logic block in 
the sixth cycle (if such a signal was formed during the given iteration), and az(tj) = —1 signals are emitted during 
the seventh cycle, The numbers Y and R are again fed to the registers from outputs S of the summators after the 
operations have been performed. Operations with negative numbers are done with an inverse code, The signals 
ay(tj) and ay(tj) which enter the integrator during one iteration are stored in the input block, The structure param- 
eters F, = NM and F, = 1 correspond to a digital analog computer consisting of such integrators, 

SUMMARY 

This paper describes a method for choosing the optimum structure of a specialized computer based on the 
derivation of a figure of merit. In. one or another way, the choice of this figure is always based on intuitive 
considerations, and so the method presented cannot give the best of all possible systems in the strict sense of the 
word, However, it gives an approximation to the best structure, the degree of which depends on a number of 
factors, and to a large degree on experience in designing such systems, 
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STABILITY OF MULTISTAGE MAGNETIC AMPLIFIERS 
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Translated from Avtomatika i Telemekhanika, Vol. 22, No. 1, 
pp. 97-106, January, 1961 
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A multistage magnetic amplifier is considered here as a linear system with time lag, The fre- 

quency criterion is applied in investigating the stability of amplifiers with rigid negative feedback, 
The stability conditions, expressed in terms of the circuit parameters, are used in the design of de- 
cision amplifiers with assigned accuracy, Experiments showed that the dynamic idealization of the 
magnetic amplifier was correct, and that the simplifications made in the theoretical analysis were 


permissible, 


Multistage magnetic amplifiers with rigid negative feedback often become unstable. In this, the system 
passes into the region of nonlinear static characteristics, and self-oscillations arise in the system, The main cause 
of instability naturally cannot always be determined, For certain given values of the time constant and the feed- 
back factor, the cause of instability is the lag in individual amplifier stages. 


Until recently, in considering the dynamics of ordinary (slow) magnetic amplifiers with time constants 
much larger than the feed-voltage period, the lag caused by the discrete character of changes in the output volt- 
age of magnetic amplifiers with internal feedback was usually neglected, This simplification can lead to sub- 
stantial errors, For instance, if lag is neglected, a two-stage 
magnetic amplifier with rigid negative feedback must always 
be stable, However, it often happens that the operation of such 
amplifiers is unstable, 


The present paper is concerned with the determination of 
stability conditions and the clarification of related problems in 


(Uo, 


\ 


if the design‘of multistage amplifiers which have magnetic ampli- 
/ fication stages with internal feedback, It is assumed that each 
/ stage operates on the linear portion of its static characteristic. 
a f We shall neglect the influence of the coupling of stages through 
r the feed source, and consider that the feed source is sufficiently 
Fig. 1. powerful, 


Dynamic Characteristics of Single-Stage Magnetic Amplifiers 


Magnetic amplifiers with internal feedback (self-saturation) are characterized by the separation in time of 
the process of core induction variation caused by the signal action — the control half-period— from the core “in- 
terrogation™ process — the operating half-period, During the operating half-period, a voltage corresponding to 
the core magnetization level in the preceding control half-period is produced at the amplifier output. In other 
words, the output voltage lags behind the input signal by one half-period of the feed voltage, 


The output voltage can be approximately characterized by averaging its actual values for each feed-volt- 
age half-period T/ 2; 


\ Up (t)dt (n=0, 1, 2,3,...)- 
nT /2 


In this, the amplifier transfer characteristic will be given by a step curve (Fig. 1, Curve 2), The envelope 
(Curve 1) of the step curve represents the transfer characteristic of the amplifier control circuit. 


N 
Yo The amplifier transfer function which would take 
| = 10s Ne into account the discrete character of the output voltage 


as | | variation can be obtained by writing the amplifier differ- 
LAI ence equation, and by applying the discrete Laplace trans- 
Fig, 2. formation to this equation [1,2]. However, in considering 


magnetic amplifiers with deep negative feedback, we can 

use Only an approximate representation of the transfer func- 
tion by utilizing the fact that, after feedback is provided, the amplifiers become basically fast-acting, It is well 
known that fast amplification stages have a lag equal to the feed-voltage half-period, Therefore, the dynamic 
properties of a stage can be given by the approximate transfer function 


Uo (p) K 


where T ; is the control circuit time constant, which, under steady-state conditions, is related to the voltage gain 
Kyj of this circuit by the well-known relation® [1] 


Kore (3) 


Here, 1 ; is the amplifier load circuit efficiency, and w.j and w..j are the numbers of turns in the control 
and the operating coils, respectively; Ktjj = Uo/ Ug, is determined by taking into account the voltage drop across 
the signal source internal resistance; f is the feed frequency. 


If tT; >> T for individual stages, the step transfer charactert:tic of a stage can be approximated by the trans- 
fer characteristic of the amplifier control circuit, which is shifted by T/4 (Fig. 1, Curve 3), Then the approxi- 
mate transfer function will be given by 


(2) 


Kee PT /4 
= 


The maximum possible lag in the load circuit was used in Eq. (2), Therefore, the approximation used can 
result in a certain increase in the stability margin, 
Method of Investigating the Stability of Multistage Magnetic Amplifiers 


The transfer function of an N-stage magnetic amplifier without feedback is determined by the expression 


N 
—pT 1 
W (p) = Kuve P ‘Wt (4) 


N N 
where T, = and Ky = || 
j=1 j=1 


Here, Tj and T; are the lag and the time constant, respectively, and Ky; is the voltage gain of the _jth stage, 


If the amplifier is encompassed by rigid negative feedback (Fig. 2), with the feedback factor 6 (for instance, 

with respect to voltage), the amplifier transfer function assumes the following form: 
Ww K 
BK y (1 + 

The necessary and sufficient conditior for the stability of a system with feedback is that the characteristic 

equation of the system must not have roots w.th positive real parts: 


* Equation (3) can be rigorously proved for nonreversive magnetic amplifiers; it is valid for the majority of re- 
versive circuits, However, if, in the case of any complex circuit, there are doubts whether it can be applied to a 


reversive amplifier as a whole, it can be applied to each single-cycle element separately. 
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1 + BW (p) = 0. (6) 


For determining the stability conditions for closed-loop lag systems, we shall use the method based on the 
application of frequency criteria [3,4]. 


Assuming that p = jw in (6), we shall obtain the equation for determining the critical frequency and lag 
time values corresponding to the stability boundary: 


BW (jo) = —1, 
or, according to (4), 
+ 
j=1 
whence we find 
BK 
N =f (8) 
Visor 
i=1 
and 
+ (m=0, +1, +2,...), (9) 
or 
T; (m=0, +1, +2....), 
1 
= (11) 
HVi+o 


j=1 
After performing calculations by using (10) and (11) for different values of w, we can plot the stability 
boundary in the form given by BKy = f(T). 
If the T; and the Kj, values are assigned for a certain given amplifier, the critical value 6¢, of the feed- 
back factor corresponding to the stability boundary can be determined from the graph obtained in this manner. 
In this case, the amplifier stability condition will be expressed by the following inequality: 


B< Ber- (12) 


If the T; and 6 values are known, the critical gain Ky, of an open-loop amplifier can be determined in 
a similar way, and the stability condition can be expressed in the following form: 


Ku < Kocr. (13) 


Finally, if the 6 and Ky values are given, the adduced equations determine the critical lag T; ._, and the 
amplifier will be stable for 


Th (14) 


Any of the above stability conditions can be used equally in calculations, As will be shown below, they 
all make it possible to relate the stability of magnetic amplifiers to their structural parameters, 


Expressions (10)-(14) are valid for amplifiers having any number of stages, as well as for mixed amplifiers, 
where not all the stages are magnetic, In the latter case, in determining the value of T; in Eq. (11) and (14), 
it is obviously necessary to take into account the fact that, for a nonmagnetic stage, the lag can be different 
from T/2, If it is necessary to take into account the effect of lag T, in the magnetic-amplifier control circuit 
(see [1], Section 6B), it is sufficient to add this value to the over-all amplifier lag. 


Ky Ky 


of 


of simplified formuja (19) 
a 
0 20 50 
Fig. 3. 


As an example, the stability of two-stage magnetic amplifiers will be considered in more detail below, 
This example is of considerable practical interest because such amplifiers are widely used, and also because, 


without taking into account the lag, it is impossible to explain the case of unstable operation of these amplifiers, 
which is observed in practice. 


Investigation of the Operating Stability of Two-Stage 


Magnetic Amplifiers with Negative Feedback 


We shall first find the critical feedback factor from Eqs, (10) and(11), which, in the case of a two-stage 
amplifier, are written as follows: 


V (1 + (1 + = BKu, Ku, (15) 


oT; = x(2m + 1)—arctgot,—arctgwt, (m=0, +1, #2,...). (16) 


The latter equation can be readily reduced to the following form: 


@ (T+ 


oT = arc tg (17) 


By solving Eq, (15) with respect to w*, and considering that the amplification factors are large, we shall ob- 
tain the approximate value 
2 _ Ky 


Wor = (18) 
which, after substitution in (17), makes it possible to obtain the following simplified expression for the critical 
value of the feedback factor: 

U; U; 1 
or, if we assume the two-stage amplifier lag time to be T; = 1/ f, 


%) 
Ky, Ky," 


The simplified equation (19) provides a convenient relationship between the amplifier lag, the time con- 
stants of individual stages, their gains, and the negative feedback value, 


In order to determine the limiting values of the ratioa® =T,/T, of the time constants of the stages and 
the limiting values of the Ky,Ky,8 product, for which the application of the simplified equation (19) secures a 
theoretical accuracy of 10% in comparison with calculations where (15) and (17) are used, we performed calcula~ 
tions for the range of & values from 1 to 50, and for the range of Ky Ky,§ vaiues from 1 to 2000, The boundary 
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curve corresponding to a 10% accuracy in calculations with respect to (19) divides the diagram in Fig. 3 into two 
portions; the lower portion of the diagram corresponds to the totality of @ and Ky,Ky,6 values for which calcu- 
lations should be performed by using (15) and (17). As we see from the diagram, formula (19) secures sufficient 
accuracy for the majority of cases encountered in practice, 


If the amplifier consists of finished stages with known voltage gains Ky, and Ky,, and the time constants 
T , and T, determined according to Eq. (3), the actual stability condition for the given amplifier can be deter- 
mined according to (12) and (19) if 8 is expressed in terms of the feedback circuit parameters. For instance, 
for an amplifier with magnetic feedback, 


Ra Wiw 
B Riw Wo,’ (20) 


where Rg, is the resistance of the first-stage control circuit, Rfy is the feedback circuit resistance, and Wfy and 
Wey are the number of turns in the feedback windings and in the input stage control winding, respectively. There- 
fore, we obtain the following stability condition from (12) and (19) for this case: 


Mag (21) 
Rew > (Ti + Ta) / Re, Wo, 


By taking (3) into account, this condition can be written in a different form: 


(22) 


Expression (22) is convenient for determining the magnetic amplifier stability by using the static character- 
istics of individual amplifier stages. 
Application of Stability Criteria in the Design of Decision Amplifiers 


Assume it is necessary to design a decision amplifier with the voltage amplification factor K, where the 
relative variation of this factor y = dK/K, i.e., the amplifier error, does not exceed ¥9. The feed source fre- 
quency and voltage , the signal-source resistance, and the load resistance are given, 


We shall assume that, from the experience gained in utilization, or as a result of corresponding measure- 
ments, the possible variations of the amplification factors of individual stages for certain given operating condi- 
tions (a given range of changes in feed voltage, frequency, and temperature) are known, i,e., that 7, = dKy, /Ky, 
and 7, = dKy, /Ky, are given. The values of y, and y, are independent of the control circuit parameters, 


The initial relationships to be used in designing are obviously given by the inequalities 


B <Ber- (23) 
The amplifier error is related to the errors of the stages by the equation [5] 
tit 


BKy, Ky, 


As a rule, for a decision amplifier, y << 7, + 72 must hold, which is secured if 8Ky,Ky, <<1, In this case, 


~ 
If we take (24) and (19) into account, conditions (23) assume the following form: 
at U+t 
Ko, < To» B< Ko, Ky, . (25) 


By combining these conditions into a single dual inequality, we obtain 


11+ %2 
BKu, Ky, < (26) 
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2Ky, Ky, Re We 
Rtyw > | | 


re Ung Condition (26) determines the range of BKy,Ky, values 
Ro for which the designed amplifiers will satisfy the necessary ac- 
curacy requirements while remaining stable at the same time. 


K Wiw We shall express (26) in terms of the assigned gain value 
K for an amplifier with feedback, According to (5), we can 
| | write the following expression for the gain: 


K Ky, Ky, 
T+BRy, Ky,’ 


or, approximately, K * 1/8, Hence, we obtain for (26): 
K 11+ Ku. K 
uv, Ku, << (27) 


Considering that 


d= 


we can obtain the following inequality from (27): 


K i 

< Ko, <K (1+ =): (28) 

by using this inequality we can determine the allowable range 

Fig. 4 of gain values for one stage with respect to certain chosen param- 
fie eters of the other stage (in this case, the second stage), 


It can happen that the necessary requirements cannot be satisfied by the available amplifier types [in this, 
inequalities (23) become incompatible], In this case, it is necessary either to make the requirements less severe 
with respect to accuracy and/or with respect to the decision amplifier gain, or to provide special measures for 
improving the amplifier operating stability, Among such expedients are increasing the feed-source frequency, 
increasing the time constants of the stages (for instance, by means of shorted windings), the introduction of sta- 
bilizing (flexible) feedback, etc, It should be noted that improvements in the quality of magnetic materials and 
valves would make it easier to satisfy conditions( 26), which would facilitate the design of stable and highly ac- 
curate decision magnetic amplifiers capable of securing considerable gain values, 


Experimental Investigation of the Stability of 


Two-Stage Magnetic Amplifiers 


For checking the obtained theoretical results, we used two two-stage reversive amplifier circuits with in- 
ternal feedback in each stage; the circuits are shown in Figs, 4 and 5, In the first circuit (Fig. 4), which operates 
at a feed frequency of 500 cps, both stages are based on a differential circuit with the median point in the load 
circuit, The displacement in the first stage was provided by shunting the feedback rectifier by means of resistor 
Reh, and separate windings, to which a suitable displacement current was fed, were provided in the second stage, 


In the second circuit (Fig. 5), which operates at a feed frequency of 50 cps, the first stage is based on a 
circuit with ballast resistors, and the second stage is based on a differential circuit with a single common load 
and commutating resistors in the operating circuits [6]. 


We tested eight magnetic amplifiers based on the above-described circuits; the amplifiers had different in- 
put and operating circuit parameters (see table), For feedback in the first stage of each amplifier, a winding 
with Wfy = 100 turns was provided, 


The voltage gains were determined with respect to the static characteristics of the amplifiers, and the time 
constants of the stages were calculated by using Eq. (3), It should be noted that the time constants must be cal- 
culated by taking into account the presence of all additional circuits, which cause an increase in the time con- 
stant, as well as by taking into account all the resistances in the control circuit, including the signal-data 


~ Ry | 
2 We, Rr 
3 R5 
4 2 
|| 
Fig. 5. 
: First-stage data 
Feed 
500 4 450 100 19.8 48 | 0.39 39 =| (0.034 
500 2 450 100 65 48 0.39 11.8 0.015 
500 -. 450 200 39.6 48 | 0.39 30.8 0.053 
500 4 450 200 85.2 48 | 0.39 11.6 0.041 
50 5 450 100 19.8 48 | 0.39 30 
50 6 450 100 1520 48 | 0.39 0.024 | 0.009 
50 7 500 1500 223 | 600 | 0.47 100 30 
50 8 500 1500 500 | 600 | 0.17 45 13.5 
Feed Secondestage data ~ Rew, ohm 
frequency, calcu- | meas- 
cps "ee Ros Ku, lated | ured 
500 1500 200 250 0.53 | 173 | 0.047 5600 5200 
500 1500 200 250 0.53 | 173 | 0.017 8500 8600 
500 1500 200 250 0.53 | 208 | 0.0163 | 3700 3200 
500 . 4500 200 250 0.53 | 208 | 0.0163 | 4450 .4300 
50 850 500 250 0.53 | 625 | 0.1468 3900 3800 
50 850 500 250 0.5° | 625 | 0.168 4250 4100 
50 300 600 50 0.53 83 6 310 270 
50 300 600 50 0.53 83 6 575 510 


transmitter resistance, For instance, for the circuit shown in Fig. 4, in the first stage of which the displacement 
is secured by shunting the internal feedback rectifiers by means of resistor Rgh, the value of T , obtained for the 


Cy s 
must take into account the internal active resistance of the first-stage output circuit, which reduces the second- 
stage time constant 7 , times, 


given stage by using (3) must be increased (1 + + Fe) mes For the second stages of all amplifiers, one 


Then, the theoretical values of the critical resistance Rew cr Of the negative feedback circuit in each ampli- 


fier were calculated in correspondence with (22), In calculations, it was assumed that the lag of each stage was 
T/2. 


The critical value of Rey, Was experimentally determined for the case where no signal was present at the 
amplifier output, The resistance Rg, in the negative feedback circuit was gradually reduced until self-oscilla- 
tions arose, This resistance value was taken as the critical value, and was compared with the theoretical value 


determined from (22), The amplifier data and the results obtained by calculation and measurement are given in 
the table, 


The experimental results indicate that the simplifications made in the theoretical analysis of the problem 
[the approximate representation of the transfer function (1) and the simplifications in deriving (19) and (26)) 
are permissible, Almost in all cases the theoretical values of Rey ¢, exceed the measured values, which confirms 
the assumption that the simplifications made lead to a certain improvement in stability, 


The satisfactory agreement between the theoretical and the experimental results indicates that, for active 
loads, the lag caused by the discrete character of changes in the output voltage of magnetic amplifiers with in- 
ternal feedback is the main cause of the instability of two-stage amplifiers of this type. 


SUMMARY 


1, In stability investigations, a multistage magnetic amplifier which consists of stages with self-saturation 
and is provided with negative feedback can be considered as a linear system with lagging feedback, The lag of 
each stage can be taken with sufficient accuracy to be equal to the feed-voltage half-period, 


2, By using the frequency criterion in magnetic amplifier stability analysis, it is possible to obtain a simple 
relation between the stability condition and the circuit parameters, and to use these conditions in rational design, 


3. For active loads, the lag of magnetic amplifiers with self-saturation is the main reason for the excita- 
tion of two-stage amplifiers with feedback, 
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SELECTING THE PARAMETER CORRELATIONS ' OF TWO TYPES 
OF THIRD-ORDER SINGLE-LOOP AUTOMATIC CONTROL SYSTEMS 
WITH ADDITIONAL PULSES WITH RESPECT TO THE DERIVATIVE 
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Translated from Avtomatika i Telemekhanika, Vol, 22, No. 1, 


pp. 107-110, January, 1961 


Original article submitted April 20, 1960 


Recommendations concerning the selection of parameter correlations of two types of third-order 
single-loop automatic control systems are given in this article, 


1, Statement of the Problem 


Third-order single-loop automatic control systems (ACS) with additional pulses with respect to the deriva- 
tive of the controlled variable are often encountered in practice, In the case of unit step excitation, the transi- 


arctgs 


Fig. 1. 


Kap 


ent processes in such systems are determined by the values of four independ- 
ent parameters while, in most cases, it is desirable that the value of one of 
these parameters — the system over-all gain— be larger than those of the 
other three, Investigations performed by several authors [1,2] indicate that, 
for third-order ACS and unit step excitation, the optimum quality of the 
transient process is secured if the ACS characteristic equation has triple real 
roots, 


2. Graphical Method for Determining _ 


Aperiodic Stability Conditions 


The characteristic equation of a single-loop ACS consisting of three. 
aperiodic elements with an additional pulse with respect to the derivative 


has the following form: 
(p + 1) (Tip +1) (Tap + 1) + Bp +X =0, (1) 


where p = d/dr; T, and T;, are the relative time constants of elements in 

the control circuit (reduced to the time constant of the system to be controlled), 
and T is the relative time (reduced to the time constant of the system to be 
controlled), 


The characteristic equation (1) can be represented by the following system of equations: 


D; (p) = (p + 1) (Tip + 1) (Tap + 1), 
D, (p) = Da (p). 


The graphic solution of the system (2) is given in Fig. 1, which, for a certain given value of K, shows the 
limiting straight lines OA and OB determining the boundaries of aperiodic stability, and the straight line CD de- 
termining the maximum value of Kap for aperiodic stability. 7 

It can be readily seen that the values of K,, and 6 shown in Fig. 1 correspond to the boundary of aperiodic 
stability pertaining to the case where the characteristic equation has triple real roots, 


Thus, the use of this graphic method makes it possible to find the values of 8 and Kap with respect to 
certain given T, and T, values, and also to find the regions of K and 8 values corresponding to aperiodic stabil- 


ity. 
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3. Conditions for the Extremum of K,, 


The straight Line CD is tangent to the D,(p) curve at the inflection point, which makes it possible to find 
the analytical expressions for 6 and Kap by using the well-known mathematical relations for third-power poly- 
nomials [3): 


ap 27 (3) 


The dependences of Kap and 6 on T, and T, are given in Figs. 2 and 3, An investigation of expression (3) 
for the extremum indicates that, if T, and T, are variable, Kap has the minimum value Kgn® for Ty = T; = 1. 
If T; = const and T, is variable, the K,,, value is obtained if the relationship between T, and T, has the follow- 
ing form: 


27's 
(5) 


The values of K3p" for the case (5), as well as the corresponding values of 8 and T, in dependence on T;, 
are shown in Fig. 4. 


4. Selection of Parameter Relations in the Case Where the System to Be 


Controlled Is Not Self-Equalizing 
In this case, the application of a similar method yields the following expressions: 


Kap = (6) 
B= (1) 
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The minimum values of K,,, are obtained for T, = var and T, = var if Ty = T, = 1; for T, = const and T, = 
= var, the minimum values of Kap are obtained if 


= (8) 


5. General Regularities for Third-Order ACS with Derivative Input 


The above method can be applied to third-order ACS with input with respect to the derivative if the ACS 
characteristic equation is given by 


ap’ + bp?+cp+d+Bp+ K =0. (9) 
The critical gain (corresponding to the stability limit) of such an ACS is determined by the expression 
b 


where Kg is the critical gain of the original ACS (without input with respect to the derivative), 
The conditions to be fulfilled if the characteristic equation (9) is to have triple real roots are given by 
B =ba—c, (12) 
where & = b/3a is the absolute value of the triple real root. 


For the value of B obtained from (12), the following relation holds: 


Kor 


SUMMARY 


1, In order to obtain triple real roots of the characteristic equations of the considered ACS types, it is 
necessary to determine the Kap and 8 values by using Eqs. (3), (4), (6), and (7), 


2, From the various obtained relationships between ACS parameters, one should exclude the relationships 
between T, and T, that are close to those given by Eqs, (5) and (8), which will make it possible to improve the 
static accuracy of ACS for the optimum quality of the transient process, 


3. The necessary (but not sufficient) condition to be fulfilled if the ACS characteristic equation of the 


type given by (9) is to have triple real roots is the condition that the ratio of the gain to the critical gain must 
be equal to 9 (for a closed-loop system), 
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A NEW TYPE OF PNEUMATIC COMPUTERS, I 


G. T. Berezovets, V. N. Dmitriev, and A. A, Tal* 
Translated from Avtomatika i Telemekhanika, Vol, 22, No, 1, 

pp. 111-118, January, 1961 

Original article submitted July 24, 1960 


This article is concerned with methods for the design of pneumatic analog computers that are 
based on pneumatic decision amplifiers operating at low pressures (0-100 mm of water), The 
first part of this article describes methods for performing linear algebraic operations, 


1, Introduction 


Earlier, the development of industrial pneumatic automation followed only along the lines of designing 
various regulating devices and servomechanisms, In recent years, communications on the development of 
pneumatic analog and digital computers have been appearing more and more frequently in our country and 
abroad [1-6 and others}, These simple and reliable devices can be successfully used in the rather large fields of 
technology where these devices are not required to operate at high speed, 


At the present time, there are many different types of components and devices, which is due to the absence 
of a unique method for designing pneumatic computers, 


There is also another difficulty to be overcome in the design of accurate pneumatic computers, namely, 
the fact that the throttle discharge characteristics are nonlinear for the standard range of input and output pres- 
sures, which is equal to 0-1 atm excess, 


Both these difficulties can be eliminated by using two new principles in the design of pneumatic devices, 


The first of these principles requires that all control and computer devices in pneumatic automation be de- 
signed on the basis of a single standard unit — the decision (operating) amplifier — as is done in the case of elec- 
tronic computers, 

This design method makes it possible to work out a unique method for designing different types of pneum- 
atic devices, and to minimize the number of necessary parts, 


| 
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Fig. 1. 


The second of these two new principles requires that the range of operating pressures be only one-hun- 
dredth of the range presently used, i.e., that pressures in the 0-100 mm of water range be used.* 


The switch to low operating pressures does not entail a deterioration of the dynamic characteristics of these 
devices; it provides a solution to the problem of designing linear throttles, and the air consumption is thereby re- 
duced at least tenfold, 


The main purpose of the present article is to demonstrate how most diverse computers can be designed by 
using pneumatic decision amplifiers operating in the 0-100 mm of water pressure range, The present article is 
concerned with different components for this system and methods for performing linear algebraic operations, The 
methods for the design of integrating and differentiating computers, as well as computers for some nonlinear 
operations, will be given in the second part of the article, 


Striving to encompass the greatest number of different computer operations, the authors will here consider 
only layout problems related to the operating principles of computers without dealing either with technical pro- 
blems connected with the construction of the devices, or problems involving their operating accuracy. 


2. The System Components 


According to the proposed principles for the design of pneumatic computers of the type under consideration, 
standard parts are used in the layout of each computer, The various parts used in the system can be divided into 
two groups, The first group pertains to the basic unit in the system, the pneumatic decision amplifier, and the 
second group comprises pneumatic throttles and pneumatic capacitors ordinarily used in pneumatic automatic 
devices, Moreover, individual units used in pneumatic relay technology (which are also standardized) (see [5}), 
can be applied ix: devices of this type. 


The pneumatic decision amplifier, whose layout and external view are shown in Fig. 1, consists of a frame 1 
with a fixed nozzle n; a membrane unit 2, with two identical thin membranes m, the rigid centers of which are 
connected by means of a bracket f , which acts as a flapper in front of the nozzle; and a permanent throttle 3 
which can be pulled out, The membrane unit is placed inside frame 1, whereby three chambers a, b, and c are 
formed, Pressures P, and P}, are fed to chambers a and b, respectively, and chamber c, where the nozzle and the 
flapper are mounted, is open to atmosphere, One side of the pull-out throttle is connected to the feed line, where 


*The idea of designing devices operating at low pressures belongs to Ferner (East Germany), On the basis of this 
principle , he constructed devices which yielded satisfactory results [1], 
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the pressure is equal to Ps..4, and the other side is connected to nozzle n. For the amplifier, P, and P), are in- 
put pressures, and the pressure P, which is established in the pull-out throttle chamber, provides the amplifier 
output signal, It is obvious that P = Ps..qg when the flapper completely covers the nozzle, and that P = 0 atm ex 
if the nozzle is not completely covered by the flapper, This is due to the fact that, instead of an ordinary, a 
pull-out nonvariable , throttle is used in this amplifier. 


It can be readily seen that the following function can be simulated in this decision amplifier in correspond- 
ence with its layout: 


P 
P 
Sg (P, — Po) (1) 


ies, P = Pheeg if Pa > Pp and P = 0 if Pz < Py (Sg is a function equal to 1 if P, — Py, > 0; it is equal to 0 if 
P,— Pp < 0), The static characteristic shown in Fig, 1c illustrates this dependence, The greater the accuracy 
with which function (1) is simulated in the decision amplifiers, the greater the accuracy of devices based on these 
amplifiers, 


In other layouts, the decision amplifier will be conventionally represented in the form shown in Fig, 1d. 
Three types of pneumatic throttles are used in this system: 


1) a nonvariable linear (laminar) throttle, whose static discharge factor is given by G = Ap, where G is 
the discharge by weight through the throttle, Ap is the pressure drop in the throttle, and & is a constant proportion- 


e greater the accuracy of ie. on the 
a (turbulent) throttle, ‘hone factor is given by G 


8) a variable linear throttle, whose discharge factor is the same as that of the comesponding senvariable 
throttle, and which differs from the latter only by the variability of the a coefficient, 


Thr 


3. Simulation of Linear Algebraic Functions — 


p, where G &s 
meaang. for, simulating a dependence of the following form: 


bi c (tur actor siven = ; 


throttle. w e tec hare: th that of the corresponding nonvariabile 


whees are the input preumatic pressure), ond P is the output 
“Tet us Consider, in tir, three methods by means Of which at problem can be solved in verlous degrees,° 


eral 


ans « Figure,2 shows, the, schematic, layout of a nosoompenseting throttle commotion device, 
In ‘this layout, the input pressures P;, Pz, «++» Pp are fed to nonvariable linear throttles &,,%2,...,%n- The 


output pressure P is established in the line eonnecting the outputs of all throttles, 


(9) 


*Each of these methods corresponds to an electronic method which has been applied to pneumatics, 
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If we neglect the output line capacity, we obtain 


>) & =, (3) 
i=1 


where G; = o;(P;— P) is the discharge by weight through the ith throttle, and o; is the proportionality factor of 
this throttle, 


From Eq, (3) we obtain 


n 
P=) kP, 
i=1 


n 
O<k <1, (5) 
i=1 


Conditions (5) considerably restrict the capability of the summation device under consideration: it cannot 
be used for ordinary addition, since ky < 1; neither can it be used for subtraction, since kj > 0. 


However, summation devices of this type, being inordinately simple, are widely used in pneumatic auto- 
mation devices for the solution of problems satisfying conditions (5), As an example, we can indicate the sum- 
mation device of the positive feedback chamber in the ACS control unit [7]. 


- ‘The Second Method, Figure 3 shows the layout of a compensating inverting throttle summation device, 
which is based on the application of a decision amplifier. 


The layout consists of n + 1 linear throttles (%,, %:,..., %p and &* ) with connected outputs and a single 
decision amplifier. Input pressures P;, Pz, Pn are fed to n throttles %2,...,%n), and pressure P, 
which arises at the decision amplifier output, is fed to the input of the (n + 1)th throttle («* ), The constant pres- 
sure P, = P® = 50 mm of water is maintained in the chamber a of the decision amplifier, and the chamber b is 
connected by means of the line connecting the outputs of all the n + 1 throttles, If we take into account the 
fact that the presence of the negative feedback shunting the decision amplifier secures the maintenance of the 
pressure Pp = P* in the b chamber, and if we neglect the capacities of the connecting line and of the decision 
amplifier b chamber, we obtain 


n 
a* (P— P*) + Dia, = 0. (6) 
=1 


By introducing the notation 
i., by measuring all pressures from the “conventional zero P* ,” we obtain 


P=— > k,P,, 


i=1 


a 


where 
k= 
and, consequently, 
(8) 
where 
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Fig. 4, Layout of the device capable 
of performing any linear algebraic 
operation by means of two decision 
amplifiers, 


In this case, the value of k; varies within the interval 


0< ky < (9) 


Any linear algebraic function can be simulated by using these summation devices, In order to demonstrate 
this, let us consider the dependence (2) when written in the following form: 


P=)} a,P, — (10) 
t=—1 i=1 


If notation (7) is used, the dependence (10) assumes the following form: 
n m n m 
P= >) a,P,— >} t > a— (145 
i=1 i=] t=1 i=1 
Moreover, considering that, if any Pj = 0, Py = —P*, the dependence (11) will be written thus: 
m+1 


P = > a,P;— > Ons 


i=l 


m n 
i=1 i=1 


Relation (12) can be written in its final form: 


mri 


> a;P,, P=— (?, > > (14) 
i=1 


The layout shown in Fig. 4 can simulate the dependence (14) and, consequently, the initial function (10), 
This confirms the possibility of obtaining any linear algebraic function by means of this layout if, in the general 
case, two decision amplifiers are used, 


par 
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Fig, 8. 


, From among the important particular variants of the lay- 
oH out given in Fig. 3, and of the corresponding dependence (8), 
the following should be mentioned, 


A, An inverting amplifier (Fig, 5) is obtained from the 
general layout if only a single input P, is used in it, Assuming 
that n = 1 in (8), we obtain the following relation for this lay- 

layout of the compensating throttle out: 
summation device of the second P =— kPy, 


type. 


(15) 
where 


B, We obtain an inverter instead of the inverting amplifier if 


a 


In this case, we have 


as) 
- The Third Method, Figure 6 shows the layout of a compensating throttle summation device of the second 
type. It differs from the previous type by the fact that the chamber a of the decision amplifier is here connected 
to the output of a noncompensating throttle summation device with m inputs, By using relation (4) for determin- 

ing the pressure P, that arises in the decision amplifier chamber a, we obtain 


m 
1 
i=1 


In the same manner, for the amplifier b chamber, we have 


(er 


a’ + > a; i=1 


The presence of compensating feedback in the layout (the decision amplifier output acts on pressure Pp 
through the &* throttle) leads to the fact that P, = P}, and, consequently, 


(19) 


a A 
Hence, 


It follows from relation (20) that the summation device under consideration makes it possible to perform 
any linear algebraic operation given by (10) by using only a single decision amplifier, in contrast to the summa- 
tion device based on the previous layout, where at least two decision amplifiers had to be used for this purpose, 


In particular, it also follows from (20) that, if 


m=n+i, Oy = (é=1,..., n) and =O", (21) 


q=1, == (i=1,...,m”) and Gi, = 1, (22) 


n 
Pa Ponti + > a, 
i=] 


A 
a 
p* p* p* 
b > BS a? 
a, a; 
A p* 
| 
4 
?) 
| 
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where 
n 
i+ 
ants 
=1 
then 
and, consequently, 
where 
99 


The layout corresponding to this case is shown in Fig. 7, This variant of the compensating throttle summa- 
tion device of the second type is often used, In the case where n = 0, this variant is transformed into the layout 
of the so-called repeater (Fig. 8), which performs the following operation: 


P= Py. (24) 


Pneumatic automation devices which perform certain algebraic operations can always be based on com- 
pensating summation devices (inverting or second-type devices) and they can sometimes be constructed on the 
basis of very simple noncompensating summation devices if their limited capabilities resulting from this are ac- 
ceptable, In this, the computers do not have to be based on summation devices of the same type; it is some- 
times advisable to use different summation methods in a computer, For an illustration of this, we shall consider 
possible layouts of proportional regulators which, as is well known, are intended for the performance of the fol- 
lowing operation: 


P =k(P;— P2) + Ps, (25) 


where P is the regulator output pressure, P, and P, are the pressures proportional to the present and the assigned 
values of the controlled variable, respectively, P; is the pressure determining the so-called control point, and k 


is the regulator gain, 


Figure 9a shows the layout of a proportional regulator based on a single compensating summation device 
of the second type, It represents a particular case of the layout shown in Fig. 7, While it is extremely simple, 
this layout has a considerable disadvantage: the adjustment of the gain value is here performed by means of two 
pneumatic throttles, 


Figure 9b shows the layout of a regulator based on four inverting compensating summation devices (two of 
them are inverters), In this layout, the adjustment of k is performed by means of a single pneumatic throttle; 
however, this results in a much more complicated layout. Finally, Fig. 9c shows a regulator layout with a single 
compensating inverting summation device and a single compensating summation device of the second type. Here, 
the gain is adjusted by means of a single throttle, and, moreover, this layout is much more efficient that that 
shown in Fig. 9b from the viewpoint of the amount of equipment used, The above example indicates that, in 
certain cases, the application of combined summation devices of different types is the most efficient procedure, 
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EXPERIMENTAL DETERMINATION OF THE DYNAMIC CHARACTERISTICS 
OF PNEUMATIC PIPES 


A. M. Smirnov 

Translated from Avtomatika i Telemekhanika, Vol. 22, No. 1, 
pp. 119-120, January, 1961 

Original article submitted January 18, 1960 


This article provides a brief description of the equipment designed for the determination of the 
dynamic characteristics of pneumatic pipes by means of the frequency method, The experi- 
~ mental frequency response curves for two pneumatic pipes are given, 


Pneumatic pipes (usually metallic pipes containing compressed gas) are used in pneumatic regulators and 
control systems, For a correct estimate of the performance of regulators under dynamic operating conditions, it 
is necessary to take into account the piping dynamic characteristics, 


( 3 


Fig. 1. Schematic diagram of the test arrangement, 


An exact theoretical determination of these characteristics is a very complicated problem, Therefore, an 
experimental method for determining the frequency characteristics of pneumatic pipes is given below. 


For the experimental determination of the frequency characteristics of pneumatic pipes, it is necessary to 
have equipment for generating harmonic pressure oscillations at the pneumatic pipe input and for measuring and 
recording pressure oscillations at the pneumatic pipe input and output without significant dynamic errors, 


The test arrangement schematically shown in Fig. 1 can be used for this purpose, The stand consists of a 
special device for generating pressure oscillations and an inductive two-channel apparatus for measuring the pres- 
sure Oscillations, 


The device for generating pressure oscillations consists of the following parts: a disk 1 with profiled open- 
ings, which is driven by an electric motor 2; a compressed-air tank 3; a control pressure gauge 4, The stream of 
compressed air supplied from the tank to the input of the pneumatic pipe under investigation, 9, is interrupted by 
disk 1, due to which near-sinusoidal pressure oscillations arise in the pneumatic pipe, This device makes it pos- 
sible to generate pressure oscillations with frequencies of up to 1000 cps and an oscillation depth (double ampli- 
tude) of up to 1 atm for a pressure of 15 atm in the pneumatic feed pipe. 


101 


| | 


The rotating disk has a diameter of 18 cm 
J and is provided with 12 profiled through-cuts, 
/ The inside diameter of the pressure-feed pipe 
and of the receiving pipes is equal to 6 mm, A 
F VV specimen of pressure recording is shown in Fig. 2. 
The two-channel inductive apparatus consists of 


Fig. 2, Sample of the recording of pressure oscillations two inductive pressure-data transmitters 5 and 6, 
with a frequency of 70 cps, 1) Pneumatic pipe inlet; which are mounted at the pneumatic pipe input 
2) pneumatic pipe outlet, | and output, an amplifier-converter 7, and a loop 


oscillograph 8, This equipment makes it possible 
to measure and regulate pressure oscillations with frequencies of up to 1000 cps with a frequency response error 
not greater than 5%, 


The inductive data transmitters have a small pressure chamber, the volume of which does not exceed 0,5 
cm’, due to which the volumes of the pressure chambers of the data transmitters that are connected to the 


100 600 700 weps 


Fig. 3, Frequency responses of pneumatic pipes: 1) pipe length 100 cm; 2) pipe 
length 25 cm, 


Thus, harmonic pressure oscillations of the form p, = pcoswt are generated at the pneumatic pipe input, 
Pressure oscillations of the same frequency w arise at the pneumatic pipe output; however, they have a different 
amplitude, pA(w), and a certain phase shift y (w); they are described by p, = pA(w)cos[wt — ¢(w)]. By com- 
paring the ordinates on the pressure oscillation recordings at the pneumatic pipe input and output, the frequency 
response 


Aw) = f Ww) 
is obtained on the oscillograph strip chart, 


Figure 3 shows the amplitude responses of two 4X 6 mm copper pneumatic pipes with lengths equal to 
100 and 25 cm, which were determined experimentally according to the above method, 


From the data in Fig, 3, it is clear that the frequency responses of pneumatic pipes in the frequency range 
of up to 700 cps have a number of resonance maxima: The pneumatic pipe 100 cm long has five maxima, and 
the pneumatic pipe 25 cm long has two maxima, The largest frequency response value occurs at the first reso- 
nance maximum, where it attains 1.7-1.8, For the subsequent resonance maxima, the amplitude response values 
gradually decrease and, at frequencies close to 600 cps, the amplitude response is equal to 0,5-0.6. 


pneumatic pipe can be neglected, 
20 
Ala | | | 
y ol 
102 


THE DESIGN OF RELAY CIRCUITS WITH THERMORESISTORS 


V. F. Bakhmut-skii 

Translated from Avtomatika i Telemekhanika, Vol, 22, No, 1, 
pp. 121-124, January, 1961 

Original article submitted June 6, 1960 


A method for designing relay circuits with thermoresistors for an assigned temperature control 
range is described, Practical recommendations and formulas are given, 


1. Temperature control devices where the relay effect is used in circuits with a thermoresistor (TR) are 
widely used in industry because they are simple and reliable, 


At the present time, two basic variants of these devices have been developed: 
1) devices with small units and fixed temperature control setting (for instance, UTS-1,[1], p. 166); 
2) devices with large units with variable temperature control setting (for instance, UTK-6 [1]). 


The following question now arises: Is it possible to design a third variant which would combine the ad- 
vantages of the above two variants, i,e., which would satisfy the following requirements: 


a) The device contains standard interchangeable relay units and thermal data transmitters, 


b) The device contains an element for temperature control setting in a given range (and, correspondingly, 
a setting scale) for each group of control points corresponding to a single mechanism or to a single portion of the 
technological cycle, 


In order to solve this problem, we shall calculate a relay circuit with TR for an assigned temperature con- 
trol range. 


2, Consider the circuit with a thermoresistor R;, an additional resistor Rg, and a correcting resistor R, 
(see figure), Let 0* be a certain temperature of the medium surrounding TR, which lies within the limits of the 
assigned temperature control range 0, = @* = 62, For determining the circuit parameters corresponding to the 
advent of the relay effect at a given temperature, we shall use the well-known graphic plot [1). 


The volt-ampere characteristics of the circuit elements are shown in the figure, where Ry denotes the yolt- 
ampere characteristic of the basic TR (TR with the largest voltage of the volt-ampere characteristic maximum), 
Ry, is the volt-ampere characteristic of the TR for which Rg is selected, I, and Up (py and Up,) are the current 
in, and the voltage across, Ry (Ry) at the beginning of the relay effect, I, and Ue are the current in, and the 
voltage across, Ry at the end of the relay effect, and Uf is the circuit feed voltage, The value of le is determined 
by the relay actuation current Ie = Igy under the condition that the relay winding resistance notexceed Rg. The 
values of R,, Ug, and R, can be found graphically or by calculation according to Eqs. (123), (124), and (133) (1): 


Up (0°) — Ve (8*) 
Ra (6*) Ie = lb )» (1) 
U; (8°) = Up (0°) + Ra (6°), (2) 


If the indicated calculation method is applied to the entire 9, = @= 0, range, we shall obtain different 
values of Ry, Ug, and R, for different @ values; in this case, for switching over from one controlled temperature 
to another, it will be necessary to change three parameters simultaneously, which is not acceptable. The pro- 
blem consists in changing one parameter so that the other two remain constant, Since Re depends on Ra, the 
values of these two resistances should be kept constant, and the temperature setting within the given range should 
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be secured by varying Ug; with such an arrangement, the relay units and the thermal data transmitters will be in- 
terchangeable, i.., the requirement a in paragraph 1 will be satisfied; it is obvious that the requirement b will 
also be fulfilled. 


We shall now determine the actual values of Rg and Re, and the dependence U¢@), Considering that the 
minimum value of I, corresponds to the maximum temperature, we shall write 


Up (82) — Ue (82) 
Ra (0) = Ry (0s) = (4) 


The determination of R, is restricted by the following requirements: TR must not be heated to a tempera- 
ture exceeding the maximum allowable temperature; the long-duration current passing through the relay wind- 
ing after it has been actuated must not exceed the maximum allowable value, The first requirement is satisfied 
by choosing a suitable relay (for TR of the KMT~-11 type, IgyUay(0,) = 0.2 w must hold] and by blocking TR when 
it is actuated, The simplest way of meeting the second requirement is to provide an arrangement for disconnect- 
ing the circuit feed after actuation or all the circuits of a given group pertaining to a single controlled mechan- 
ism, In the case where it is necessary to change the setting and to continue the controlling operation after the 
circuit has been actuated, the feed is temporarily disconnected, The time during which the feed is to be dis- 
connected is determined experimentally from the condition that there should be no spurious actuation: This time 
must be sufficient for TR to cool down (approximately) to the surrounding medium temperature after it has been 
heated to 300-400°C at the moment of actuation; delay time must be effected automatically by means of a time re- 
lay. Besides disconnecting the relay, its overloading can be prevented by other means, For instance, in the case 
where the number of changes in the feed voltage U¢(@) in the 0, = @ = 0, range is close to (or less than) the num- 
ber of times the actuation and dropout currents pass through the relay, it is advisable to introduce ballast resistors, 
which are connected in parallel or in series to the relay winding when it is actuated; this would make it possible 
to reduce the relay winding current to the maximum allowable value for the largest feed voltage U(9,) while 
keeping its value larger than the relay dropout current for the minimum feed voltage U¢(@2). 


The values of Us@), where 0, = © = 6,, can be obtained graphically (if a series of volt-ampere character- 
istics of Ry is available), or analytically by substituting in (2) the value of Rg calculated from Eq, (4), and the 
1,,(@) and U},@) values found from the general volt-ampere characteristic equation [2] 


U B 
©) 
where A, B, and F are the parameters of Ry. 


With an accuracy sufficient for practical purposes, we can consider that 1,(@) and Up@) coincide with the 
coordinates of the volt-ampere characteristic maximum point; by using the condition for the existence of maxima 
dU/dI = 0, we obtain from (5) 
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Ub a hy 
Ry 
U 
|| 


8 V +) B 
Ip (0) = A exp | Up (8) 


a (0) = —0 +2-y (-8 +BY - 


For RT, Eqs, (5) and (6) are of the same form, but the A,, B;, and F, values are different, 


Let us now formulate the conditions under which R, can be kept constant throughout the entire 9; = 6 = 6, 
range. By substituting (6) in (3), it can be readily shown that Rg depends on @ (for a constant Rg) only due to a 
spread of the B values or, which is the same, due to a spread of the TR temperature coefficients, while R, in- 
creases with a decrease in@. Consequently, by assuming that Rc = Ro (02) for B, # B, a temperature error is in- 
curred, which is the larger, the higher the value of 6, in comparison with ®,;, The maximum allowable spread 
of B is found from the following condition (cf, (3)): 


(0: + Abr) 


where R¢ is the value of R, corresponding to the 
MO, > 0 is the maximum allowable error, After substituting in (7) the expression for R,(@,) from (3) and the ex- 
pressions for Uj and Up, from (6), and after performing a number of calculations where the higher powers and 


products of the small quantities “1—* ¢, and are neglected, we obtain the following value for the maxi- 
1 
mum allowable spread of B at the temperature 0;: 


B, —B Af, 

[ (8) 
If Ro = Re), the error will be of the opposite sign and, instead of (8), we obtain 

Bi — — Ads 

(9) 


where A@, <0 is the maximum allowable error for the temperature setting 0,, and [55], =the comespond- 
ing spread of B values, 


It is obvious that the error will be smaller throughout the entire range if R, remains within the range of 
Re = Rc) values, where 0, <0) < 43. Similarly to (8) and (9), we obtain, in this case, 


From (10) and (11) we obtain 


0: +0 [Bi—B) _ 200. 


The temperature coefficient spread of KMT-11 thermoresistors from a batch with the same rating usually 
does not exceed 0.3% per °C, which corresponds to a “t— spread of up to 0,06, If we assume that Ps 
= 0,06 for the range 35-100°C (@, = 308°K, 0, = 373°K), we obtain from (12).0)~ 68°C,A@, ~ 2°C, Thus, by 
using the R, value for 68°C, the setting error can be reduced to + 2°C throughout the entire 35-100°C range, with- 
out additional corrections for the spread of TR parameters; the error sign is determined during calibration by 
suitably shifting the setting points, 
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where 
Re (03) = Re (01 + AQ;) = | (7) 


It should be noted that the setting error 46 constitutes only a portion of the over-all error; its other por- 
tion consists of errors due to an inaccurate selection of R, and R,, an inaccurate adjustment of the feed voltage 
Up @), changes in the load, etc, The latter errors also arise in the case of fixed setting and, under mass-produc- 
tion conditions, they can be reduced to + 2-3°C; in this case, the over-all error will be 4-5°C, In certain cases, 
the indicated error may be unacceptable; then, in order to reduce AO y, it is necessary either to select a TR with 
respect to the temperature coefficient, or to average the TR characteristics by connecting in series several TR 
with a lower rating in a single thermal data transmitter. 


8. On the basis of what has been said, the following conclusions can be drawn; 


1) Thermal control devices, where the relay effect in circuits with TR is used, can be provided with ad- 
justable actuation temperature setting (with an error of up to +5°C in the 35-100°C range) while preserving the 
miniaturization of the component units and the interchangeability of the relay units and thermal data transmitters, 


2) The devices embodying these properties are flexible; they do not contain a great amount of wiring, and 
they are convenient in use and adjustment. 


3) The construction of these devices does not entail the use of circuits and designs which would be greatly 
different from those presently available [1]; however, the requirements pertaining to the spread of TR parameters 
or, more accurately, the spread of TR temperature coefficients, must generally be more strict, 


In conclusion, the author extends his thanks to G. K, Nechaev for the discussion of the paper, and his valu- 
able advice, 
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EVENTS 


SEMINAR-CONFERENCE ON THE THEORY AND METHODS 
OF MATHEMATICAL SIMULATION 


G. M. Kozyreva 
Translated from Aytomatika i Telemekhanika, Vol, 22, No, 1, 
pp. 125-126, January, 1961 


The first session of the Seminar on the Theory and Methods of Mathematical Simulation, which was organ- 
ized on the initiative of the Institute of Automation and Remote Control, Academy of Sciences of the USSR, was 


held on July 12-13, 1960 in Moscow, Problems pertaining to the design and application of digital simulators 
were discussed at the session, 


More than 60 people from 17 scientific-research organizations in Moscow, Leningrad, and Gor*kii partici- 
pated in the work of the Seminar, Five reports were heard and discussed, 


In his introductory speech, B. Ya, Kogan noted that the Institute of Automation and Remote Control, 
Academy of Sciences of the USSR, considers that the solution of the question of organizing the discussion of the 
most important problems in computer techniques and, especially, in the field of their application in actual auto- 
matic systems, has been long overdue, In connection with this, the Institute intends to hold annual all-Union 
seminar-conferences, It is intended that each of these conferences be devoted to a certain important problem, 
and that it summarize the results achieved during the past year and indicate the direction of further development, 
B. Ya. Kogan further indicated that a situation has now arisen where the accuracy and the dynamic range of elec- 
tronic simulating devices is insufficient for the solution of a number of important problems in the automatic con- 
trol field. These limitations of electronic simulators are especially noticeable in solving complex (high-order) 
nonlinear differential equations, in the case where it is necessary to perform multiple coordinate transformations, 
and in the simulation of functions of a single and, especially, many variables, 


The utilization of universal digital computers for these purposes is not justified because of the complexity 
of equipment and programming. 


In recent years, a new method — digital simulation — has been developed, which apparently can be of great 
help in overcoming these difficulties, 


A. V. Shileiko submitted a report entitled "A method for synthesizing the optimum structure of digital 
simulators," In his report, he proposed a method for synthesizing the optimum structure of specialized computers, 
which consists in the selection of the performance criterion, the determination of its dependence on a number of 
parameters, the so-called structural coefficients, which determine the structure of the device realizing an as- 
signed algorithm, and in finding all the structural coefficient values securing the maximum performance criterion 
value, The amount of information provided by the device in unit time, reduced to the device complexity, is 
used as the performance criterion, This method was used for synthesizing the optimum structure of an actual 
digital simulator, 


The report by L, M. Goldenberg, entitled “Sequential digital differential analyzers (DDA)," was devoted 
to the description of the DDA “Integral and provided some of its technical characteristics, The machine is of 
the sequential type and has a ferrite core memory. It has 32 integrators, The cyclic-pulse frequency is 50 ke, 
The iteration duration and, consequently, the integration speed, depend on the number of integrators fitted and, 
if their number is complete, it is 20 iterations per second, The machine contains 350 tubes and 8000 ferrite cores, 


The report by K. S. Neslukhovskii, entitled “Integrating adapter to UTsVM,” described a sequential-action 
DDA intended as a special adapter to the BESM-2 machine for solving differential equations, A standard mag- 
netic drum memory is used in this device, The adapter operates at a rate of 50 iterations per second, In order 
to speed up the DDA operation, the programming of initial conditions is performed by means of a universal digi- 
tal.computer, 


In his report, “Sequential DDA for higher frequencies," F. V. Maiorov reported the development of sequen- 
tial DDA for operation in control systems, 


In his communication, entitled "Frequency characteristics of DDA," S. V. Misailovskii proposed a method 
for determining the transmission band of DDA in dependence on their parameters, and he also determined the 


maximum actual rate of change in input signals produced without dynamic errors for a given frequency trans- 
mission band. 


The following resolution was adopted at the seminar: 


1, The initiative of IAT AN SSSR in organizing a permanently acting seminar on the theory and methods 
of mathematical simulation has been approved, 


2. It is considered necessary to intensify efforts in designing digital simulators, to systematize the ex- 
perience gained in their utilization, and to specify the fields of their possible application, 


3, It is noted that problems pertaining to the methods and theory of digital simulation are not treated to 
a sufficient extent in the domestic technical literature, 
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